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Chapter I 



Introduction 



A description of a large system of particles is often sought in a derivation from 
the detailed behaviour of just a few of the particles. The present thesis deals with 
the connection between such microscopic features and the nature of a collection 
of many particles. A study of identical bosons is an obvious first investigation of 
this hnk, but the ideas might also be apphed to fermion systems or to systems 
with mixed symmetry. The relatively well-known, flexible properties of cold 
alkali gases have presented questions which might be addressed by a study of 
few-body correlations. In this chapter wc comment on basic features of bosonic 
systems and motivate a description of few-body correlations within a many- 
particle system.* 

1.1 Bosons 

All particles can be classified as either bosons or fermions. The distinction is 
important when identical particles approach each other. Electrons, nucleons, 
and quarks are fermions and obey the Fermi-Dirac statistics, while force car- 
riers like photons and gluons are bosons and obey the Bose-Einstein statistics. 
Atomic nuclei, atoms, and molecules obey one of the statistics depending on 
the number of contained fermions. In this thesis we as far as possible consider 
bosons generally, but often relate to bosonic neutral atoms with an even number 
of neutrons, in particular alkali atoms like ^''Rb or ^"^Na. Although bosons are 
the main objects, we will a couple of times consider extensions of the methods 
to deal with fermions. 

At most one fermion can occupy the same quantum state, whereas bosons 
are not restricted. An example is Bose-Einstein condensation of a vast number 
of bosons in the same single-particle state, which was experimentally achieved 
in 1995 by coohng dilute alkah gases [AEM+95, BSTH95, DMA+95]. 

A mean field is the basis for the description of dilute alkali gases by the Gross- 
Pitaevskii equation [EB95, BP96]. In mean-field self-consistent theories [BJ83] 

*The use of "we" refers to the author with the reader's participation. 
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the influence from interactions is included as an average, hence the name mean 
field. Such a description is reasonable when the interaction between particles is 
so weak that each particle only feels the other particles as an average background 
cloud in which they move. The rigorous criterion is that the mean free path is 
long compared to the spatial extension of the system. Reviews of theoretical 
developments before and after the experimental realization of Bose-Einstein 
condensation arc given in the references [DGPS99, PS02, PSOSj. 

At strong interaction or large densities each particle might interact strongly 
with one or a few of the other particles. The particles then adapt to the local 
environment, but still feel the (weaker) mean-field influence from the remaining 
particles. This competition between the background and the local surroundings 
is important when we formulate the theory in chapter 2. 

When the attraction is too large, for example when condensates collapse 
[SSH98, SGWH99, DCC+01, RCC+01], correlations beyond the mean field are 
crucial. Descriptions based on the mean-field Gross-Pitaevskii equation [Adh02b, 
US03] can account for this collapse, but avoid direct relations to underlying mi- 
croscopic processes. 

Some methods go beyond the mean field, but still avoid the explicit inclusion 
of correlations. An example is the Skyrme-Hartrcc-Fock method with density- 
dependent interactions [SJ87]. Related are low-density expansions of the total 
energy for a many-boson system [LY57, BHM02]. 

For large densities particle encounters are more frequent and at least two- 
body correlations need to be explicitly included. This cannot be done directly on 
top of the usual mean field with a two-body contact interaction, i.e. of zero range, 
which would lead to a wave function with zero separation and diverging enorgj^ 
[FJOla]. On the other hand, the use of realistic potentials in self-consistent 
mean-field calculations leads to disastrous results because the Hilbert space does 
not include correlations as needed to describe both the short- and long-range 
asymptotic behaviour [EG99]. 

An explicit inclusion of correlations is done by the Jastrow method [Jas55] 
where the many-body wave function is written as a product of two-body am- 
plitudes instead of one-body amplitudes. With a few assumptions about the 
asymptotic behaviour of the amplitudes, this results in variational numerical 
procedures that can be carried through for many-boson systems also for large 
densities [MMOl, CHM+02]. 

1.2 Two-body properties 

A study of the properties of a many-particle system requires understanding of 
the two-body problem. The interaction between neutral atoms is repulsive at 
short distances and attractive at large distances. There may for alkali atoms 
be a large number of bound two-body states that are sensitive to the details 
of the interaction. However, when two atoms approach each other slowly from 
afar, their encounter can be described in a universal way irrespective of the 
short-distance details of the interaction. The determining parameter is then the 
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scattering length. This is infinitely large in the presence of a two-body bound 
state with vanishing energy. Then the two particles correlate in all space, which 
is the opposite limit than assumed by a mean field where no correlations are 
allowed. 

At a Feshbach resonance, when the energy in a scattering channel coincides 
with the energy of a bound state in another channel [PS02], the scattering 
length diverges in the same way as when a two-body bound state occurs. In 
recent years such a Feshbach resonance has been investigated for both sodium 
[IAS+98, SIA+99] and rubidium gases [RCB+98, RBC+01]. For ^^Rb atoms 
an external magnetic field can slightly change the effective two-body potential 
curves, which have resulted in a tool for tuning the two-body scattering length 
[CCR+00]. 

Mean-field studies of dilute boson systems often include a two-body contact 
interaction with a coupling strength given by the scattering length. At large 
densities this becomes a problem since the interaction energy for the bosons 
then diverges. It is also difficult to handle large scattering lengths close to a 
Feshbach resonance. The alternative use of a boundary condition given by the 
scattering length at zero separation between the bosons allows larger density 
and scattering length [FJOlbj. 

Ultimately, the best description of the physics properties could be obtained 
by using realistic potentials which besides the correct large-distance two-body 
behaviour also incorporate high-energy features and the correct nature of bound 
states. However, this is especially difficult when the two-body system contains 
innumerable bound states, as is the case for the alkali atoms in experiments. 
Furthermore, if the goal is a description of the low-energy two-body properties 
within the many-body system, it would be an investment of too much effort 
in the wrong place. A more rewarding method is to use a simpler finite-range 
potential with the correct scattering length, for instance a linear combination 
of Gaussian potentials [BGOl]. Considerations about the two-body interaction 
return in chapters 3 and 5. 

1.3 Few-body physics 

A system with only a few particles can be described accurately without crude as- 
sumptions. The related methods can provide insight into complicated problems, 
for instance how two or three particles approach each other within a many-body 
system. If the spatial extension of the system is large, an encounter of two partic- 
les can be considered as a pure two-body process with an average background 
influence from the other particles. Alternatively, it might be considered as a 
three-body process, the third body being the collection of the other particles. 
In smaller systems three particles approach each other more frequently, which 
then demands a description of a true three-body process. Faddeev [Fad60] 
wrote a wave function as a sum of terms that account for pairwise encounters, 
which by Yakubovskii [Yak67] was extended to account for the right behaviour 
of three-body clusters, four-body clusters, and so on. Within each cluster it is 
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then possible to treat the degrees of freedom rather accurately, while keeping in 
mind that the present collection of particles moves relative to the particles out- 
side the cluster. The art of these Faddeev-Yakubovskii techniques is to make the 
proper assumption about the dominating structure of the many-body system; 
otherwise little is accomphshed and the calculations turn out as comphcated as 
for some methods based on a Jastrow ansatz. 

Specifically for bosons, but implemented in nuclear physics, an approach 
with inclusion of two-body correlations was worked out by de la Ripelle et 
al. [dlR84, dlRFS88]. This approach is equivalent to the Faddeev-like equations 
to be described in chapter 2. The result is an eigenvalue equation in only one 
variable. Barnea [Bar99b] proposed the inclusion of higher-order correlations in 
a method which reminded of the Yakubovskii technique. 

The hypcrspherical adiabatic method, which was formulated for a study of 
the helium atom by Macek [Mac68] , separates the description of the three-body 
system into a common length scale, i.e. the hyperradius, and an additional 
hyperangle. The extension of this method is now frequently used in atomic 
physics for descriptions of many-electron systems [Lin95] . 

After the experimental realization of Bose-Einstein condensation, most de- 
scriptions of this phenomenon were based on the mean-field Gross-Pitaevskii 
equation. However, Bohn et al. [BEG98] introduced a hyperspherical method 
to study the effect of interactions within the many-boson system. This hyper- 
spherical method was simplified by the inclusion of a contact interaction and the 
assumption that the only dependence is on the average distance from the centre 
of mass, i.e. no correlations were allowed, which reminds of a mean field. An 
advantage of the hyperspherical method is that it provides an effective potential 
in a linear eigenvalue equation, in contrast to the non-linear nature of the Gross- 
Pitaevskii equation. In a related study Blume and Greene [BG02] calculated the 
properties of three bosons in an external trap with no assumptions about the 
structure of the wave function, and thus confirmed some of the behaviours of 
the effective potentials in the hyperspherical model for the many-boson system. 

The detailed study of three particles provides an important first step in 
the formulation of a theory for clusterizations within a many-body system. 
The Faddeev-formalism is often applied within the hyperspherical adiabatic 
approximation when the emphasis is on the asymptotic two-body properties 
[JGF97, NFJGOl]. The threshold phenomenon of infinitely many bound three- 
body states in the case of a two-body bound state with zero energy [Efi70] can be 
described by just the inclusion of two-body correlations [FJ93]. This anticipates 
that a generalization of the method to a many-boson system might describe the 
case of large scattering length where non-correlated models become inadequate. 

1.4 The thesis work 

The work for this thesis started from a few-body description of two-body cor- 
relations within a many-body system, and it has been centred on solving the 
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many-boson problem in a hyperspherical frame. Central questions are formu- 
lated as follows. 

What is the effect of two-body correlations? 

The outset for the thesis work was to understand how two-body correlations 
influence the properties of a many-boson system. This work is closely related to 
studies of the three-body system, see e.g. Nielsen et al. [NFJGOl], and to the 
hyperspherical investigation of the average properties of the many-boson system 
by Bolm et, al. [BEG98]. The formulation of the main techniques behind the 
inclusion of such correlations was given in the publications [SFJN02, SFJ02b]. 
This is collected in chapter 2 and appendix B. The effects of two-body corre- 
lations were mainly discussed in the references [SFJ02b, SFJOSa] and are here 
collected in chapters 3 and 4. 

What happens close to a resonance? 

The case of large scattering length provides for the three-body case the Efimov 
phenomenon of many bound three-body states [Efi70, FJ93]. The possibility of 
similar threshold effects for the many-boson system was investigated during the 
thesis work and the results were pubhshed in [SFJ02a]. This is discussed mainly 
in chapters 3 and 4. 

Are the deviations from the mean field trustworthy? 

The relations to the mean field and the deviations of the results were published 
in [SFJ04], where also ranges of validity were considered. This is here included 
in chapter 2 and in chapter 5. 

Do two-body correlations influence stability? 

The macroscopic stability problems for a many-boson system were briefly dis- 
cussed in the previous publications, but were further addressed in relation to 
the observed phenomenon of macroscopic collapse [RCC+01, DCC+01] in the 
reference [SFJOSa], which also included a discussion of the competition between 
three-body recombinations, macroscopic tunneling, and macroscopic collapse. 
Chapter 6 collects these considerations. 

How does the geometry influence the system? 

Deformation effects were discussed in a preprint [SFJ03b], which gave stability 
criteria, effective dimensions, and effective interactions for bosons in a deformed 
external field. This discussion is continued in chapter 7. Since an inclusion of 
correlations in the deformed case is not presently implemented, this treatment 
is somehow similar to a mean-fleld treatment. 
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1.5 Thesis outline 

The structure of the thesis can be summarized as follows. Chapter 2 presents 
the hypersphcrical frame for studying few-body correlations within the many- 
boson system. This includes a discussion of the structure of the wave function 
and some comments on the mean-field wave function. We present the necessary 
assumptions when restricting to two-body correlations and the main features 
of the resulting equations of motion. Chapter 3 contains a discussion of an- 
alytical and numerical solutions to the hyperangular part of the Hamiltonian. 
Chapter 4 contains a discussion of the hyperradial part of the description of 
a many-boson system, as well as a discussion of Bose-Einstein condensation in 
the hyperspherical model. In chapter 5 we compare results from the previous 
chapters to the results from the mean-field Gross-Pitaevskii equation and ranges 
of validity are estimated. In chapter 6 we discuss stability criteria, approaches 
to the dynamical evolution of the many-boson system, and relevant time scales. 
Chapter 7 deals with the effects of a deformed external trap. Finally, chapter 8 
contains conclusions and discussions of the results. 



Chapter II 



Hyperspherical description of correlations 



Hyperspherical methods are used in studies of both few-body and many-body 
problems, for example in atomic physics [Mac68], nuclear physics [JGF97], and 
within atomic physics especially for many-electron systems [Lin95]. Especially 
relevant in the present context is a study of a many-boson system in a hyper- 
spherical frame performed by Bohn et al. [BEG98]. A large number of degrees 
of freedom are described in terms of one length, the hyperradius, and some 
angles called hyperangles. This reminds of the characterization of a two-body 
system by the relative two-body distance and the angular degrees of freedom. 
In the hyperspherical description, the angular degrees of freedom are usually in- 
tegrated or averaged such that the many-body system is described by only the 
hyperradius. An effective potential depending on the hyperradius then carries 
the information about the average angular properties. This is once more anal- 
ogous to the description of the two-body problem by only the radial distance 
with inclusion of the angular momentum in an effective centrifugal potential. 

In section 2.1 we first define a set of coordinates appropriate for a study 
of correlations in the many-body system. The steps are similar to the ones 
written by Barnea [Bar99b, Bar99a]. For clarity we formulate this in three 
spatial dimensions, but it can as well be written in lower dimensions. We then 
rewrite the Hamiltonian and the Schrodinger equation according to this choice of 
coordinate system in section 2.2. Appendix A contains further details. As basic 
input we need an anzatz for the many-body wave function. We discuss the basis 
for assuming a Faddeev-like decomposition of the wave function and relate this 
to other common approaches in section 2.3. We also discuss possible extensions 
to fcrmion symmetry and three-body correlations. Wc arc then equipped to 
solve the angular equation in section 2.4, which is done with the inclusion of 
two-body correlations both by a Faddeev-like equation and by a variational 
equation. Some details of the derivations are given in appendix B. We end the 
chapter by considering the complications involved in an angular equation with 
the inclusion of three-body correlations. 
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2.1 Jacobi vectors and hyperspherical coordinates 

The system of N identical particles may be described by TV coordinate vectors 
and momenta p.,, labeling the particles by the index i = I, .... N. Here a more 
suitable choice of coordinates is the centre-of-mass coordinates R — XliLi fil^ i 
the N — 1 relative Jacobi vectors rjf. with 

j N-k ( 1 \ 

= V Ar3fcTTi^^-^+^ - E (2-1) 

and k = 1, 2, ... , N—1, and their associated momenta. These Jacobi coordinates 
are illustrated for the first six particles in figure 2.1. The notation is r]k = |t7/;|, 

1 
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Figure 2.1: Jacobi vectors connecting the first six particles. 

so r?jv-i is proportional to the distance between particles 1 and 2, riN-2 is 
proportional to the distance between particle 3 and the centre of mass of 1 and 
2, riN-3 is proportional to the distance between particle 4 and the centre of 
mass of the first three particles, and so on.* 

Hyperspherical coordinates are now defined in relation to the Jacobi vectors. 
One length, the hyperradius p, is defined by 

I ^ N N 

p"^PN-i = j^Y. 4 = E(^^ - ' (2-2) 

k—l i<j i—1 

where = \ri — rj\. The last two equalities show that the hyperradius can 
be interpreted either as {N — l)/2 times the root-mean-square (rms) distance 
between particles or as \/N times the rms distance between particles and the 
centre of mass. 

In three spatial dimensions, the N — 2 hyperangles ak G [0,7r/2] for k = 
2,3, . . . , N — 1 relate the length of the Jacobi vectors to the hyperradius via the 
definition 

smak = — . (2.3) 
Pk 

* Throughout the thesis normal font for a corresponding vector denotes the length of that 
vector, i.e. tj^ = \ri,,\, n = \ri\, etc. 
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Since pi = ?7i, the fixed angle ai = 7r/2 is superfluous, but is for convenience 
often included in the notation. Remaininff are the 2{N - 1) angles n'fl'^ = 
{•dk, f>k), for k = 1, 2, . . . , iV — 1, that define the directions of the N — 1 vectors 
rjj., that is 'dk G [0,7r] and tpk G [0, 27r]. All angles are collectively denoted by 
f2 = {ak, fk} with A; = 1, 2, . . . , — 1. In total the hyperangles and the 
hyperradius p amount to 3(A^ — 1) degrees of freedom and the centre-of-mass 
coordinates R amount to three. These coordinates are also related by 

N N . . N .2 

E^^ = ^E4 + ^(E-0 =p' + NR^- (2.4) 

i=l i<j ^ i=l ' 

The total volume element is 

JV N-l 

l[dri = N^/^dRl[drjk, (2.5) 

i=l fe=l 

where the part depending on relative coordinates is Ilfclr/ c^^fe-^ hypersphe- 
rical coordinates this relative part becomes 

N-l 

H drjk = dpp^""-^ dflN-i , (2.6) 

fe=i 

dnk = dn'-^'> dVl'^J;^ dilk-i , drii = df^W , (2.7) 
dO.'^^ = dak sin^ "fc cos^''"^ at , rffi^''^ = d'&k sin '&kdtpk , (2.8) 

(k) 

where dflri is the familiar angular volume element in spherical coordinates. 
Since the angle a^v-i is related directly to the two-body distance ri2 by sinaAr_i 
= r/N-i/pN-i = ri2/{V^p), the volume element in equation (2.8) related to this 
angle is especially important, that is d^i^"^^ = daN-i sin^ aN-i cos^^~'^ ajv-i- 
The angular volume integrals can be computed to 

where F is the gamma function [Spi68]. An angular matrix element of an oper- 
ator O with two arbitrary functions \1/ and $ is with equation (2.6) for fixed p 
then given by 

(*|o|*)o = j drtN-i **(p, n) 6 $(p, , (2.10) 

which in general is a function of p. 

In this section the many-body system is described by a straightforward or- 
dering of particles as {1, 2, 3, . . . , A^}. Furthermore, we use the configuration 

tThe notation is dvi = dri^driydriz with = {rix,riy,riz). 
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principle that a Jacobi vector connects one particle with the centre of mass of 
some other particles. Both the ordering of particles and the recursive config- 
urations have to be done differently when we later evaluate matrix elements. 
However, the present formulation of the basic equations is independent of such 
considerations and therefore the simplest structure is presented. 

2.2 Schrodinger equation for N identical particles 

With the preceding choice of coordinates we obtain the Hamiltonian and next 
rewrite the Schrodinger equation by an adiabatic expansion, which was first 
used in a study of a helium atom [Mac68] . 

2.2.1 Hamiltonian in hyperspherical coordinates 

We consider N identical particles of mass m interacting only through two-body 
potentials Vij = V{rjj). Wo do not consider effects due to spin, and thus omit 
explicit spin dependence throughout this thesis. 

An external trapping field Krap confines all particles to a limited region of 
space. This is written expHcitly as an isotropic harmonic-oscillator potential of 
angular frequency w, i.e. for particle i it is given by Vtrap(^i) = moj'^rf/2. This 
confining field is relevant for studying trapped atomic gases, but can later be 
omitted from the general results by putting u) = 0. The total Hamiltonian is 
here given by 



which with equation (2.4) is separable into a part only involving the centre-of- 
mass coordinates and a part only involving relative coordinates. The centre-of- 
mass Hamiltonian is 



where Pr = X^^^i Pj is the total momentum and M = Nm is the total mass of 
the system. We subtract this from the total Hamiltonian and get 




(2.11) 




(2.12) 



H = -f^total — Hr 



(2.13) 




Using equation (2.4) we can write this as 




(2.14) 
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Here T is the intrinsic kinetic-energy operator which in hyperspherical coordi- 
nates can be rewritten as 



T = 



1 



d 3^_4 d 



2m\p^^-'^dp^ dp 



A2 



(2.15) 



The dimensionless angular kinetic-energy operator is recursively defined 

by 



A2 



cos^ ak sm ak 



d'^ 3fc-6- (3fc-2)cos2afe d 



sin 2afe 



dak 



(2.16) 
(2.17) 



where Mk is the angular-momentum operator associated with r/j,. Thus, the an- 
gular kinetic-onorgy operator is a sum of derivatives with respect to the various 
hyperspherical angles. Convenient transformations to avoid first derivatives in 
equations (2.15) and (2.17) are 



2m ' 



-(3N-4) 9 r)3JV-4 9 

^ dp^ dp 



= p 



-(3JV-4)/2 



dp^ 



(37V-4)(3iV-6) 



ill = sin-iafeCos-(3'=-4)/2c,, 



4p2 



(3JV-4)/2 



g/c- 10 



+ 



(3fc-4)(3fc-6) 2 
■5^ '-i tan^ afe 



sinafecos'^^'' ^^^"^ au ■ 



(2.18) 



(2.19) 



The Hamiltonian H can now be collected as 



H 



AT 



2rV,, 



(2.20) 
(2.21) 



where Tp is the radial kinetic-energy operator, /ifj is a dimensionless angular 
Hamiltonian, and Vij is a dimensionless potential. Thus, the intrinsic Hamilton- 
ian H contains a part which only depends on p and a part hu which depends 
on n and on p through the two-body potentials Vij . 



2.2.2 Adiabatic expansion and equations of motion 

Since the total Hamiltonian is given as i^totai = Hji + H, the total wave function 
for the A'^-particle system can be written as a product of a function T depending 
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only on R and a function ^' depending on p and the 3A'' — 4 angular degrees of 
freedom collected in fl, i.e. 

*totai = T(fi)*(p,0) . (2.22) 

The centre-of-mass motion for the total mass M = Nm is determined by 

HrT{R) = ErT{R) . (2.23) 

From equation (2.12) the corresponding energy spectrum is obtained as that of 
a harmonic oscillator, that is -E_r,„ = huj{2n + 3/2) with n = 0, 1, 2, — 

The relative wave function ^{p,Q,) obeys the stationary Schrodinger equa- 
tion 

H^{p,n) = E^{p,fl) , (2.24) 

where E is the energy in the centre-of-mass system. This is solved in two steps. 
First, for a fixed value of the hyperradius p we solve the angular eigenvalue 
equation 

(/io-A,)$,(p,r!) = 0. (2.25) 

The angular eigenvalue Xv{p) depends on p. Second, the collection of angular 
eigenfunctions $i/(/9, ^l) is used as a complete set of basis functions in an expan- 
sion of the relative wave function. This is for each value of the hyperradius p 
written as 

oo 

*(p, n) = J2 Mp)Mp, fi) . Mp) = p-^^^'-^'^'^Uip) , (2.26) 

i/=0 

where the factor p^(-''^^4)/2 introduced to ehminate first derivatives in p, 
see equation (2.18). The expansion coefficients for fixed p, f^, or F,^, are then 
considered as hyperradial wave functions. 

In analogy to the technique for iV = 3 [JGF97], equation (2.26) is inserted in 
equation (2.24), equations (2.20) and (2.25) are used, and the resulting equation 
is projected onto an angular eigenfunction The result is a set of radial 
equations 



d^_2mE Kipl ^ (3jV-4)(3jV-6) V _ ^(2). ^ 
dp^ h? ^ p^ ^ 4p2 +64 ^u.KP} 



E 



'2Qi'Up)^^+Q5{p) 



U'{p) , (2.27) 



that couple the different angular channels. Here bt is the trap length given by 
bt = ■\/h/{muj), and the coupling terms Q^^li are defined as 



{<^Ap,^) 






($.(p,f^) 


MP,^))a 
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A special result is Q^J = [NFJGOl]. The angular eigenvalues Xi, enter these 
coupled equations as a part of a radial potential. The total radial potential 
Uv{p), entering on the left hand side of equation (2.27), is: 

2mU.{p) _ K{p) (37V-4)(3iV-6) 

f9 ~ p" V hi "^--^P'- ^^-^^^ 

This includes a p^-term due to the external harmonic field, a p~^ centrifugal 
barrier-term due to the transformation of the radial kinetic-energy operator, the 
angular potential A,y, and the diagonal term QuJ . 

The expansion in equation (2.26) is called the liyperspherical adiabatic ex- 
pansion. Its efficiency relies on small coupling terms Q,^^, which then requires 
inclusion of fewer channels v. In the following investigations of the dilute boson 
system, the non-diagonal terms arc often found to be smaller than 1% of the 
diagonal terms. Without these couplings the right-hand side of equation (2.27) 
vanishes, and the equation simplifies significantly to 

-|^^ + t/.(p)-i^]/.(.)=0. (2.30) 

In this thesis small couplings are generally assumed and only results of this 
non-coupled treatment are shown. 

Thus, the centre-of-mass motion is separated out and the hyperspherical 
adiabatic method turns out to be promising for a sufficiently dilute system due 
to small coupling terms. The remaining problem is the determination of the 
angular potential A from the angular eigenvalue equation. 



2.3 Wave function for identical particles 

So far no specific structures are assumed. The allowed Hilbert space for the 
many-body wave function in principle includes any structure of the system. 
However, at this point an ansatz for or approximation of the angular wave 
function $i/(p, il) is necessary. 

The Hartrec wave function with a product of single-particle amplitudes 
[BJ83] is the basis for mean-field treatments of many-particle systems. In 
contrast, Faddeev-Yakubovskii formulations [Fad60, Yak67] contain additive 
decompositions of the wave function which explicitly reflect the possible asymp- 
totic large-distance behaviours of cluster subsystems. A different starting point 
is the Jastrow factorization into products of two-body amplitudes [Jas55]. The 
Jastrow form is more efficient for large densities, while Faddeev-Yakubovskii 
methods are more successful for smaller densities where the system separates 
into smaller clusters. 

In the present hyperspherical formulation, the angular wave function can be 
written as a general expansion in the full angular space and subsequently be re- 
duced to yield a practical wave function. The result of such considerations is an 
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ansatz which is well suited for the low densities encountered for Bose-Einstein 
condensates. Two-body correlations are expected to be most important. Pos- 
sible extensions of the method to include three-body correlations or fermion 
antisymmetry are briefly discussed. 

2.3.1 Hartree: single-particle product 

The Hartree ansatz with a product of single-particle amplitudes [BJ83] is 

N 

*H(ri,r2,...,rjv) =n^H(n) ■ (2.31) 

i=l 

For the ground state of non-interacting bosons trapped by the spherically sym- 
metric external field of trap length bt , the amplitudes are given by 

VH(ri) = Ce-^^/^^'t) , = . (2.32) 

With the relation X^^Xi = + NR^ this is rewritten as 

*H(ri, r2, . . . , rjv) = exp - ^ ^ j 

= c7iVg-^V(2(.?)e-^«V(2bt=) = ro{R)Fo{p)^o , (2.33) 

which turns out as a product similar to equations (2.22) and (2.26). The separa- 
tion of the centre-of-mass motion assures that the ground-state centre-of-mass 
function always is To(i?) = CN^/'^ expl-NR"^ /{2bl)]. Then equation (2.33) is a 
product of the ground-state wave function for the motion of the centre-of-mass 
in a trap and the lowest hyperspherical wave function -Fb^o in equation (2.26), 
where -Fo(p) oc exp[— p^/(26^)] and the angular part $0 is a constant. This 
implies equivalence between a Hartree-Gaussian wave function and lack of de- 
pendence on the hyperangles ft. 

The interactions produce correlations in such a way that the hyperspherical 
wave function ^ deviates from a hyperradial Gaussian multiplied by a constant 
hyperangular part. Therefore the Hartree product wave function is strictly not 
exact. However, a measure can be obtained by calculating the single-particle 
density n from the obtained function ^, that is 




dradrs • • • (irjv|To(-R)*(p, n)f . (2.34) 



This can then be compared with the Hartree analogue |^/'h(''i)P- When the 
numerical hyperspherical solution is inserted, the 3(A^ — l)-dimensional integral 
in equation (2.34) is rather complicated. In order to get an idea of the possible 
relations, we instead assume a constant angular part $o- Then the hyperradial 
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density distribution is expanded on Gaussian amplitudes with different length 
parameters Uj 

\np)\' = E yf^N-lsN-S ^-'^^''^^ . (2-35) 

where ^ ■ Cj = 1 assures that F{p) is properly normalized as Jq°° dpp^^~^\F{p)\^ 



1. This yields the single-particle density 

3 ^ 



(2.36) 



which is equivalent to (r^) = J dry n{r\)r\ since 

kr\) = ^ip') + (R^) = l{^-^) E^.S^ + l^b' (2-37) 



and 



/ 



3 . 

dn n{ri)rl = ^^cjE"^ 



3 



The mean-square distance between the particles is then obtained by the relation 

«^>=i^(«>-(^^))=^(«>-^I^O- '^-^^ 

These relations are derived and valid only for Gaussian wave functions. How- 
ever, the true Hartree solution is not strictly a Gaussian although such an ap- 
proximation rather efficiently describes the dilute boson system [PS02]. The 
above results relate a Hartree density distribution to a similar hyperradial dis- 
tribution provided that the angular wave function is assumed to be a constant 
which corresponds to an uncorrelated structure. 



2.3.2 Faddeev-YakubovskiT: cluster expansion 

The effect of correlations is beyond a mean field where particles only feel each 
other on average and do not correlate. With the Faddeev-Yakubovskii tech- 
niques the proper asymptotic behaviours of the wave functions are directly 
taken into account [Fad60, Yak67]. These formulations are well suited when 
the large-distance asymptotics arc crucial, as expected for low-density systems. 

Faddeev [Fad60] studied three-particle systems where one of the two-body 
subsystems was bound and the other subsystems were unbound. The wave 
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function was written as $ = $12 + $13 + $23 with the three terms given by 
suitable permutations of 

$23 = </>23(r23)e'"^''^+'"^^''=^ , (2.40) 

where -R23 = (m2r2 + TO3r3)/(m2 + m^) is the centre of mass of the bound 
subsystem and Hki and ?IK23 are the momenta of particle 1 and of particle pair 
2-3, respectively. This form accounts for the details of the possibly bound pair 
2-3 and considers other effects as low-energy plane waves. A generalization of 
this three-body wave function is 

^ij = (t>ij (vij) exp M ^ KkTk + iKijRij ] , ^ = ^^ij . (2.41) 

When all relative energies are small, that is when Kij ~ and Kk — 0, the result 
is ^ij ~ <l>ij{r-,j). 

A generalization to an ^''-particle system was formulated by Yakubovskii 
[Yak67], who arranged the particles into possible groups of subsystems and 
thereby included the correct large-distance asymptotic behaviour for all clus- 
ter divisions. The decisive physical properties are related to the division into 
clusters, which for A/' = 3 amounts to three possibilities. The three Faddeev 
components are related to the number of divisions and not the number of par- 
ticles. For N > 3 the number of cluster divisions is much larger than A'^. For A'^ 
particles the wave function is therefore written as a sum over possible clusters 

*Y = X] *Y (cluster) . (2.42) 

clusters 

This method is often apphed in nuclear physics [CC98, FG02]. In a diltite 
system two close-lying particles are found more frequently than other cluster 
configurations. Then the dominating terms in the cluster expression in equa- 
tion (2.42) arc due to the two-body clusters, and the remaining particles are 
considered uncorrelated and described by plane waves or as a mean-field back- 
ground. The Yakubovskii wave function then reduces to a Faddeev-like form 
similar to equation (2.41) 

N 

vf-Y ^ Hp, ^) = Y1 ^^^^P' • (2.43) 

i<j 



2.3.3 Jastrow: two-body factorization 

The Jastrow variational formulation [Bij40, Din49, Jas55] was designed to ac- 
count for correlations in a Bose system. The Jastrow ansatz 

JV 

*j = n ^jC'^y ) ' = '^3 - '^i ' (2-44) 

i<j 



tThe complex number T is here denoted by i. 
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provides an argument for writing the wave function as a sum of two-body terms 
in the dilute limit. The two-body Jastrow component can be written as a 
Gaussian term, which corresponds to mean-field amplitudes, multiplied by a 
modification expected to be important only at small separation, i.e. 

^jK)=e-'''^/(2iV6,='[l + <^j(ri,)] , Mrij)=0 for nj > ro . (2.45) 

Beyond some length scale ro, deviations due to correlations vanish. With equa- 
tion (2.4) this leads to the relative wave function 

N 



-pV(26?) 



n[ 

i<3 



l + 0j(T'ij) 



(2.46) 



N 



N 



i<j i<j^k<l 

For a non-interacting system the sums are zero and the Gaussian mean-field 
Hartree ansatz from equation (2.33) is obtained. For a sufSciently dilute system 
it is unlikely that more than two particles simultaneously are close in space, that 
is when both nj < ro and rki < ro- Therefore the expansion in equation (2.46) 
can be truncated after the first two terms, i.e. 



N 



i<j 



~ 1 



N 



N 



N{N -l)/2 



+ 0j(^u) 



(2.47) 



A redefinition of the two-body amplitude results in a Faddeev-like sum as in 
equation (2.43). 



2.3.4 Hyperharmonic expansion of two-body components 

The sum of two-body terms can be formally obtained as the s-wave reduction 
of an expansion on a properly symmetrized complete set of basis functions. 
Appropriate are the hyperspherical harmonics y that are eigenfunctions of the 
grand angular kinetic-energy operator A^_-^, equation (2.16) [Smi60]. These 
are for collective angular momentum Lfe and projection Mfc given by [Bar99a] 



where qk denotes the set of quantum numbers {qk} = {h, 
and the hyperangular momentum is given by 



j{<]k} 



(2.48) 



Kk = 2z//s -I- Kk-i + Ik , 
The expression for y is 



(2.49) 



'{Qk} 

[Kk,Lk,Mk] 

k 



'Yiki'^k^Vk) 



n sin'^ Q,-cos^-i aj p[^+i/2.^.-i+(3j-5)/2](cos2Q,-) \ . (2.50) 



i=2 
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Here the y;,„i's are the usual spherical harmonics, Vi, is the Jacobi function, and 
the coupling of angular momenta is given by the Clebsch-Gordan coefficients 
through 



Yi, , (fii) (g) Yi^ (i?2, ¥'2) O • • • O Yi^ {dk,fk) 



^ {limil2m2\L2M2) {L2M2hm3\L3M3) •••x 



mi,m.2,---,Tnk 



(2.51) 



Omitting dependence on ak, we also define a reduced function y by 



v{9fc} 



Yi, (g)Yi^{'&2,<p2)(E>---(E> Yi^ {-dk, Vk) 



[Kk-i,Lk,Mk] 

I J[ sin'^- aj cos^^-i aj pfe+i/2,if.-i+(3j-5)/2] ^^^^ 2a,) | , (2.52) 

where {qk} = {h, . . . ,lk,V2, ■ ■ ■ i^k-i}- 

The wave function for each fixed value of p is for fixed relative angular 
momentum L = Ljv-i and projection M = M/^_i decomposed as 



N 



(2.53) 



where the component <l>y is focused on the particle pair i-j. Each of these 
components, for instance $12, can be written as the complete expansion 



{gN-i} 



[K,L,M] 



(2.54) 



where the sum runs over ah possible quantum numbers Qn-i- Analogies are 
given in [NFJGOl, Bar99a].^ The number K = Km -2 is used for reference to 
the kinetic-energy eigenvalue. No assumptions are made yet. 

In order to obtain an explicitly symmetric boson wave function in equa- 
tion (2.53), we need to rewrite equation (2.54) as 



{9JV-1} 



(2.55) 



where the second sum accounts for all possible permutations P12 of particles 
apart from the pair 1-2. 



§See page 384 in [NFJGOl] and page 1137 in [Bar99a]. 
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Zero angular momenta L = M = is reasonable in the short-range limit. 
Vanishing hyperangular quantum number K = for the remaining degrees 
of freedom yields only zero quantum numbers {gjv-i} = {0}, and the sum is 
truncated to include only the term 

^12 [0,0] (P, ^) = 4l\o,o] iP^ ^i2)i^[o°o,o] • (2-56) 

The remaining terms ^ijjti2 are obtained in similar ways, so the angular wave 
function is 

N 

^p,n) = J2<PiMnj) , (2.57) 

i<j 

omitting the constant 3^[o,o,o] and superfluous indices. 

2.3.5 Fermion antisymmetry 

For identical fermions the total wave function has to be antisymmetric under 
permutation of any two particles. If the spin wave function is symmetric, then 
the spatial wave function can be antisymmetrized by an extension of the method 
written for bosons 

N 

{«iV-l} 

The second sum of the first line accounts for all different permutations Py of 
the N — 2 particles apart from i and j. We assume that a specific permutation 
resulted in the particles ordered as ijkl.... The number p is then the total 
number of permutations of two particles needed to transform the straightforward 

ordering 1234 . . . into the ordering ijkl 

If no term depends on more than A'' — 2 particles' positions, it is not possible 
to write a properly symmetrized wave function for a system of identical fermions 
as a sum of terms. Therefore, the dependences can not be truncated as roughly 
as was done for bosons in the steps leading to equation (2.56). However, the 
antisymmetric function with lowest possible quantum numbers might provide a 
useful fermion wave function which can be implemented in calculations. In this 
thesis we do not discuss fermions further, but restrict ourselves to the case of 
bosons. 

2.3.6 Two-boson direction-independent correlations 

As seen in the preceding sections, a relative wave function of the form 

N 

^{p,n) = F{p)J2^ij{p,Q) (2.60) 

i<j 
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incorporates both mean-field properties through F{p) and correlations beyond 
the mean field through the Faddeev-components $. We therefore decompose 
the angular wave function equation (2.26), into the symmetric expression of 
Faddeev components for fixed values of the hyperradius p 

N 

<P{p,n) = J2^ij{p,n) , (2.61) 

i<j 

where each term $y is a function of p and all angular coordinates Ct. Since 
each term in itself is sufficient when all O. degrees of freedom are allowed, this 
decomposition is exact. At first this ansatz seems clumsy by introducing an 
overcomplete basis. However, the indices i and j indicate a special emphasis on 
the particle pair i-j. The component $y is expected to carry the information 
associated with binary correlations of this particular pair. 

This wave function is a natural choice for the trivial case of N = 2. A 
wave function rewritten as a sum of terms has also been successful in three- 
body computations. The advantage is that the correct boundary conditions 
are simpler to incorporate, as expressed in the original formulation by Faddeev 
[Fad60] intended for scattering. Still, mathematically nothing is gained or lost in 
this Faddeev-type of decomposition. For weakly bound and spatially extended 
three-body systems, s-waves in each of the Faddeev components are sufficient 
to describe the system [NFJGOl]. This is exceedingly pronounced for large 
scattering lengths where the Efimov states appear [FJ93, JGF97, NFJGOl]. 

The present A''-body problem is in general more comphcated. However, 
for dilute systems essential similarities remain, i.e. the relative motion of two 
particles that on average are far from each other is most likely dominated by 
s-wave contributions. Each particle cannot detect any directional preference 
arising from higher partial waves. Implementation of these ideas in the present 
context implies that each amplitude for a fixed p only should depend on the 
distance rij between the two particles. Thus, we assume 

$y {p, O) ~ (pij {p, Uij ) , (2.62) 

where the two-index parameter aij is defined by 

sin a,,- = ^ . (2.63) 

These aij's are distinctively different from the a^s of equation (2.3). 

The boson symmetry implies that the functions (pij are non-distinguishable, 
so the indices are omitted. The resulting angular wave function is 

JV N 

$(p, ^) = Y1 ^(P^ ) = Yl > (2-64) 

i<j i<j 

where (pij{p,aij) = <p{p,aij) = 4>{oLij) with the notational convenient omission 
of the coordinate p. The wave function in equation (2.64) is symmetric with 
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respect to the interchange of two particles, i ^ j, since a^- = aji and since 
terms hke 4'{aik) + 4>{ajk) always appear symmetrically. 

This ansatz of only s-waves dramatically simplifies the angular wave func- 
tion. The original overcomplete Hilbert space is reduced such that some angular 
wave functions can not be expressed in this remaining basis. Thus, the reduc- 
tion resulted in an incomplete basis, but the degrees of freedom remaining in 
equation (2.64) are expected to be those needed to describe the features of a 
dilute system. 

In section 2.3.4 the Faddeev ansatz equation (2.43) was formally established 

as a generalized partial wave expansion in terms of the hypersplierical harmonic 
kinetic-energy eigenfunctions. The two-body s-wave simplification then appears 
as a truncation of this expansion, which also leads to equation (2.64). 

In conclusion, when the system is dilute, the Faddeev ansatz with two-body 
amplitudes is expected to account sufficiently for the correlations and at the 
same time keep the mean-field-like information about motion relative to the 
remaining particles. 

2.3.7 Three-body correlations 

An extension of the inclusion of pairwise correlations to study three-body cor- 
relations in denser systems is possible and could yield insight into the process of 
three-body recombination within iV-boson systems. Since all degrees of freedom 
are kept in every term, the Faddeev-like decomposition of the wave function can 
describe all kinds of clusterizations in a particle system. The indices ij just re- 
fer to the correct asymptotic behaviour of the two-body scattering properties 
between particle i and j in a given amplitude A higher-order correlated 
wave function is in that sense included in the general expansion. However, for 
actual apphcation it does not provide any solvable method. Therefore, we pro- 
ceed as fohows with what might be an applicable three-body expansion of the 
many-body wave function. 

A symmetric boson wave function with three-body amplitudes is 

N 

i<j kjti,j 

where depends on the distances between the three particles i, j, and k. In 
hyperspherical coordinates the dependence, exemplified for the term (/'i2,3, can 
be reduced to be on p, a^-i (= ctu), OiN-2 (= "12,3), and (= ^^12,3), 

where '&N-2 is the angle between ry^_i and t;;v-2- ^ general term (f)ij^k depends 
on p, aij, cxij^k, and 'dij^k- It is written as 

4>ij,k{p,otij,aij^k,'(^i3,k) = 4>2{p,aij) + (f)z{p,aij,aij^k,'(^ij,k) , (2.66) 

where the term cix,, accounts for an n-body correlation. This is analogous to 
a proposal by Barnea [Bar99b]. Since the functional dependence is the same 
for all terms, the symmetry is explicitly included. The two-body correlated 
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method reappears with (f>s = since then (pij^k 4>2 ^ 4>{P;0/.ij). However, 4>3 
provides the three-body correlation on top of the two-body correlation. Terms 
with k < j are responsible for an overdetermined expansion since three terms, 
in principle, describe the same three-body correlation. In the case i < j < k 
the three similar terms are <pij^k, 4'ikj, and 4'jk,i- With only one term (pij^k 
we would have to impose other symmetry restrictions on this single term such 
that 4>ij,k = (i'ikj = 4>jk,i- With the sum 4>ij.k + 4>ik.j + 4>jk,i in the ansatz for 
the wave function, the symmetry is explicitly built in and is independent of the 
amplitude's functional form. A simpler description is obtained when neglecting 
■dij,k in (63, and thus yielding a term that accounts for one particle's relations 
to the pair of particles. 

2.4 Angular eigenvalue equation for two-boson correlations 

Since the eigenvalue A from equation (2.25) carries information about the two- 
body interactions and kinetic energy due to internal structure and as well about 
possible correlations, the techniques and approximations used to find A are 
especially important. 

This section contains the essential rewritings of the angular equation with 
the ansatz from equation (2.64) for two-body correlations. We first present the 
Faddeev-like equations and next construct a variational equation as an alter- 
native which is solvable under the additional assumption of short-range inter- 
actions, i.e. the system must be relatively dilute. The Faddeev-like equations 
were previously written in this form for a boson system by de la Ripelle et 
al. [dlRFS88], whereas the angular variational equation according to the au- 
thor's knowledge is an original contribution by the author and co-workers and 
first presented in [SFJN02]. end the section by briefly considering the in- 
clusion of higher-order correlations. 

2.4.1 Faddeev-like equation 

Insertion of the ansatz for the boson wave function in equation (2.64) along with 
equation (2.21) into equation (2.25) yields 



N 



N 




(2.67) 



k<l 



i<j 



with (f)ij = (f){aij). Rearrangement of summations leads to 




N r 



N 



(2.68) 



For three particles and with the assumption that each term in the square brack- 
ets separately is zero, the Faddeev equations are obtained. They have been 
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applied for the three-body case to describe particularly the regime of large scat- 
tering length [JGF97]. The same assumption for the iV-particle system results 
in the N{N — l)/2 Faddeev-like equations 



N 

(a^_i - X^(pki + Vki ^ (pij = , 



(2.69) 



i<j 



which are identical due to symmetry. A description of the many-boson system 
with such Faddeev-like equations was previously performed by de la Ripelle et 
al. [dlRFS88], who concentrated on systems within the realms of nuclear physics. 
With k — 1 and / ~ 2 the kinetic-energy operator A^_-^ from equation (2.16) 

2 



reduces to 11^ 



because A^_ 



and I 



N-l9l2 



0. Since T]j^_i = 

(r2 — r\)l\f2 and /Ojv-i = />, then equations (2.3) and (2.63) yield Q!Ar-i = O-vi- 
Therefore, only derivatives with respect to 0:12 remain, and it is convenient to 
introduce the notation II? 



12 



In the sum over angular wave function components in equation (2.69), only 
three different types of terms appear. When fc = 1 and I = 2, these types are 
classified by the set {i, j} either having two, one, or zero numbers coinciding 
with the set {1,2}. Then equation (2.69) is rewritten as 







n? 



u(ai2) 



N 

i=3 



Aj(/>(ai2) + 



(2.70) 



N 
J=3 



N 

E' 

3<i<j 



v{aki) 



V{V2p sin aki) 



(2.71) 



for a central potential V{r). Multiphcation of equation (2.70) from the left by 
(f){Qi2) followed by integration over all angular space except ai2 results in an 
integro-differential equation in a = ai2 of the form^ 

= [Ul2 + v{a)-\](t){a)+v{a)2{N-2) j dT(t){ai3) 

+v{a)^{N - 2) (AT - 3) y dr </.(a34) • (2.72) 

Here dr oc <if2jv-2 is the angular volume element, excluding the a dependence, 
with the normalization f dr = 1. Due to symmetry between the first and 
second sums in equation (2.70), this projection leaves for every value of a only 
two different integrals. Both can analytically be reduced to one-dimensional 
integrals. The results are cohected in appendix B.2. For brevity here the terms 
are denoted by 



dr (l){a34) = R^^" 



2) 



<^(a) 



(2.73) 



^Throughout we interchange a, qjv— i, and ai2. 
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I 



(2.74) 



where is an operator acting on the function (p resulting in a function of 

a. II Equation (2.72) can now be written as 



nfa + v{a) -X + 2{N- 2)v{a)RY; 



(JV-2) 



+ -(7V-2)(iY-3)«(a)R34 



(N-2) 



(2.75) 



which is Unear in the function (p. An advantage of this equation is that it does 
not become more compHcated as the number N of particles increases. 

In this equation, the potential v and the kinetic-energy operator 11^ are di- 
agonal in the sense that they only act in the space of particle pair 1-2, whereas 
the angular wave function is also evaluated for other two-body pairs. In the 
Faddeev-like equation (2.75) for the many-body case, all wave-function compo- 
nents are projected onto s waves in the 1-2 system. This means that effective 
contributions from higher partial waves in the hyperangular space are omitted. 

Another problem is that the Faddeev approximation is not variational, i.e. the 
energy may be underestimated [NFJGOl]. This is expHcitly obvious when we 
add and subtract a constant vq from the interaction potential to rewrite equa- 
tion (2.67) as 



JV N 



(2.76) 



k<l 



v'ki =Vki-vo, A' = A - -N{N - l)t;o • (2.77) 

The Faddeev approximation then results in an angular eigenvalue A which de- 
pends on the choice of Vo . This shows that the Faddeev-Hke equation has to be 
handled with care and at worst that it is inconsistent with the present assump- 
tion of s waves. 



2.4.2 Variational angular equation 

Proceeding with the Faddeev-like equation (2.75) is one option, but as discussed 
this equation shows inadequacies under the required assumptions about the 
many-body wave function. In this section we therefore rely on the full angular 
equation and discuss a variational equation where the Faddeev approximation is 
not necessary. Thus, we might effectively catch influences due to higher partial 
waves from the different subsystems in the many-body system and benefit from 
the maintained validity of the variational principle. However, in some of the 
model calculations in the following chapters, we study results obtained from the 



II Mathematically R resembles a rotation operator, hence the choice of notation. 
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Faddeev-like equation and then compare to results from the following variational 
equation. 

First the optimal angular equation is derived within the Hilbert space defined 
by the form of the angular wave function in equation (2.64). The very short 
range of the two-body interaction compared with the size of the system simplifies 
the problem as we shah see in the next section. 

The angular Sclirodingor equation for fixed p in equation (2.25) and the 
ansatz for the wave function in equation (2.64) allow the eigenvalue expressed 
as an expectation value, i.e. 



A = 



($|feo|$)n 









/n 



(2.78) 



For an operator O which is invariant when interchanging any two particles, the 
terms ((^j'j/|0|$)n = ((jf)j//j/'|0|$)o are identical since the possible differences 
vanish when averaging over all angles fi. Since ft-o is invariant with respect to 
interchange of particles, this identity holds for both numerator and denominator, 
so equation (2.78) simplifies to 



(<^i2|/in|Ei<j'/'ij>n 



yi2 



(2.79) 



The total angular volume element is (KIm-i = dVla' ^^dflif^ see 
equation (2.7). Since the integrands are independent of then dCt^~^^ 

can be omitted from the integrations. Using equation (2.10) we then obtain 



12 



dn 



N-2 



N 
i<j 







(2.80) 



The wave-function component (j)l2 is varied until the lowest eigenvalue is ob- 
tained. This gives the integro-differential equation 



JV 

/ dflN-2'^ 



k<l L 



N 



(A^_i - X)(f)ki +Wfei^< 







(2.81) 



where the unknown functions (t>jj = ct>{aij) all are the same identical function of 
the different coordinates a^ . Many terms arc identical, e.g. J dilj^^2 ^^12034 
i dQ.isi-2 Vi2<p56, since particles 1 and 2 cannot distinguish between other pairs 
of particles, see appendix B.3.1 for the details. Collecting ah terms yields 



/ 



dQ 



N-2 



(n?2 + vi2 - A) (Ai2 + G{t, aia) 



= 0, 



(2.82) 



where r denotes angular coordinates apart from ai2. The kernel G contains all 
non-diagonal parts involving other particles than 1 and 2. This is given by 



G{T,ai2) = -jn2 



1134 + t;(Q!i2) + v{aii) - A (?i(a34) 
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+ ^n2v{a34)^{oii2) + 2niw(ai3) [^{au) + 0(Q!23)] 

+2ni n^g + w(q;i2) + u(ai3) - A ^(aia) 

+n3|w(a34) [0(035) + '/'(ais)] + v {a 13)41 (a 45)^ 
+2n2v{ai3) [(/)(q;i4) + 0(a24) + (/'(a34)] 

+2n2u(a34)(/)(Qi3) + ^n4t^(a!34)?5'(Q;56) , (2.83) 

where rii = 11^=1 (-^ ~ i ~ 1) cind nf^- is defined from equation (2.17) with k = 
N — 1 and with ak replaced by a^. In equation (2.83) all terms depend at most 
on coordinates of the six particles 1-6. The first three terms in equation (2.82) do 
not depend on the integration variables r leaving only G(t, 0:12) for integration. 

By appropriate choices of Jacobi systems [SS77], the relevant dogroos of 
freedom can be expressed in terms of the five vectors tJn-i, ■ ■ ■ ,»7iv-5- One is 
the argument of the variational function and not an integration variable. The 
remaining twelve-dimensional integral is then evaluated with the corresponding 
volume element dr oc Y[i=2 d^a^ ^^dQ^^ where the normalization is J dr = 
1. Then equation (2.82) becomes 

n?2 + ^^("12) - a] (A(ai2) + j dr G{t, 012) = , (2.84) 

where the first terms are independent of the integration variables. Equation (2.84) 
is a linear integro-differential equation in one variable containing up to five- 
dimensional integrals, see appendix B.3.2. As is the case for the Faddeev-Hke 
equation, this equation does not complicate further at large A'', i.e. when 
increases beyond N = 6, the structure does not change. 



2.4.3 Short-range approximation 

The two-body potentials V{rij) are assumed to be characterized by a length 
scale b beyond which the interaction vanishes, that is when ^ b. The angular 
eigenvalue equation (2.84) simplifies in the limit when this two-body interaction 
range b is much smaller than p. Then the integrals arc cither analytical or reduce 
to one-dimensional integrals. This reduction could in principle be accounted for 
by substituting the interaction potential with a 6 function, but is done generally 
for any finite-range interaction as long as the range is small compared to the 
hyperradius. However, this is only possible for the potentials appearing under 
the integrals. Thus, apart from the local terms containing v{a), the results 
mainly depend on a parameter ae related to the volume average of the potential 
by the definition 



Tn f 



(2.85) 
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A finite value of this volume integral is essential for the validity of the method. 
This is obeyed for short-range potentials that fall off faster than . 

As an example of the reductions, when pcosa ^ b, the J dr ?;(Q!34)-term 
reduces to 



/ 



dr v{a34) = vi{a) 



2 r( 



3jV-6^ 



Ob 



TT T{^^) pcos^a 



Similarly the J dr u(ai3)-term reduces to 



J dr V 



(ttis) = V2ia) 



3\/3 



cos-'^-^/Jo e(a < 7r/3) vi{a) 



(2.86) 



(2.87) 



in the limit when p cos a cos /?o S> b, where sin/3o = tana/\/3. Here 6 is 
the truth function, i.e. it equals unity when the argument is true and zero 
otherwise. The remaining terms can in this limit be expressed through vi{a), 
V2{a), R-*^^ from equations (2.73) and (2.74), and other related operators ^[j'li- 
Corresponding definitions are given in appendix B.3.3. 

The reductions can be understood qualitatively via figure 2.2 which shows 
the geometry when the short-range interaction contributes to the integrals. In 
the integral / dr v{ai3)<p{a34), see figure 2.2a, the dominant contributions occur 
when particles 1 and 3 are close together as shown in figure 2.2b. Then the 
distance between particles 3 and 4 appearing in ^34 is approximately equal 
to the distance between particles 1 and 4. Therefore / dr v{ai3)4>{a34) ~ 
fdr v{ai3)(l){ai4)- 





034 ~ (pi4 

vi3 ~ (5(ri3) 

a) 3b)' 

Figure 2.2: Simplifications due to short-range potentials. 

So, for b <^ p the exact short-range shapes of the potential are not important 
and the integral in equation (2.84) of equation (2.83) can be written as 

dr G(t, a) ~ -^v\{a) + 4niV2{oi) (t>{oi) 



+ 



{Ri^-')n^,</,(a) + [v{a) - A]Rir'V(a)} 
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+2m{R[^-^^nl,cf>{a) + [v{a) - X\R[^-^U{a)} 
+^t;i(a)R^^"^V(a) + ^^3^^l(Q^)R^V35^(a) + n3W2(a)R^3^45 <?!>(«) 

+2n2V2 (a) [2R^^l^^(a) + R^IL '/'(")] • (2-88) 

The variational equation with these reductions is the basis for the calculations 
in chapter 3. 



2.4.4 Variational equation for three-body correlations 

The ansatz for the angular wave function from equation (2.66) includes a ge- 
neral correlation within all three-body subsystems. We choose a trial wave 
function and write the angular potential A as an expectation value analogous to 
equation (2.78) 

. JV . N N JV 

i' <j' k'jti' ,j' i"<j" i<j kjtij 

The calculation of these expectation values requires at most twelve degrees of 
freedom which with a short-range potential for 6 <C p reduces to at most nine 
degrees of freedom. 

Performing the variation 3 ~* 4^*2 s + ^0*2 s leads to the angular varia- 
tional integro-differential equation in ajv-i, oiN-2, and '&N-2- 



I 



df(A%_^ + Vi>j, - Aj ^ ^ (/>y,fc = , (2.90) 

i'<j' ' i<j kjtij 



where df denotes the angular volume element for all angles apart from ajv-i, 
aN-2, and ^n-2- There are 126 different ^-terms (38 for = 4), 12 differ- 
ent A^-terms, and 12 different A-terms. In the short-range limit many terms 
are identical and thus reduce the complications. The integrals in the integro- 
differential equation are three dimensions lower than those in the expectation 
value since three angles are fixed. Thus, the short-range approximation results 
in an integro-differential equation in three variables with up to six-dimensional 
integrals. This is beyond the scope of the present work, but indicates the com- 
plications when including higher-order correlations. 



Chapter III 



Interactions and the hyperangular spectrum 



In the hyperspherical formulation of the many-body problem in chapter 2, the 
tedious problems are "hidden" in the angular equation. The angular solutions 
carry essential information about interactions between the particles and about 
internal kinetic energy. The correlations were assumed to be two-body for sufS- 
ciently dilute systems, and this was built into the wave function. The key quan- 
tity is then the function A, equation (2.25), which determines the properties of 
the radial potential, equation (2.29). The angular wave functions potentiahy 
carry information about couplings between the different adiabatic channels. 

First, in section 3.1 we discuss how to model two-body interactions in the A''- 
particle system. Analytical derivations of angular potentials in various regimes 
are given in section 3.2. Then we comment on the numerical procedure before 
solving the angular variational equation. Section 3.3 presents the attributes of 
the found wave functions and angular potentials for various kinds of interaction 
strengths. Section 3.4 summarizes the nature of the angular potentials, which 
can be parametrized by the interaction parameters and the number of particles. 
The details behind this parametrization were previously published [SFJ03a] and 
hence collected in appendix D. 

3.1 Interactions between neutral bosons 

The effective two-body interactions vary enormously for different boson systems 
depending both on the nature of the bosons in question and on the surroundings. 
Here we consider bosons with short-range interactions in the sense that the vol- 
ume integral of the two-body interaction potential is finite. Neutral atoms, that 
are frequently encountered in experiments with dilute boson systems, interact 
via a potential of sufficiently short range and can be considered by this method. 

The interaction between atoms is repulsive at short distances due to the Pauh 
exclusion principle which forbids overlapping centres. Neutral atoms attract 
each other at longer distances due to mutual polarization which induces a dipole 
moment. The interaction between two particles, e.g. 1 and 2, can be modelled 
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by the two-body potential [GBOl] 



Kdw(r-) 



-c{r—ro) 



r = r2 - ri 



(3.1) 



This potential has the van der Waals (vdW) tail —Cq/t^ when r ^ ro, and is 
thus of short range in the sense that it decays faster than 1/r^. The important 
part of this potential is illustrated in figure 3.1 for some choice of the parameters 
Ce, To, and c. 



-10 



_ 1 1 1 1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 _ 

Geltman 

Gaussian - 









6 8 10 12 14 16 18 

r (arbitrary units) 



Figure 3.1: Two-body potentials. The solid line is the potential in equation (3.1) 
from Goltman and Bambini [GBOl], and the dotted line is the Gaussian potential 
from equation (3.7). 



At large particle separations, a direction-independent behaviour is expected, 
which means that zero relative angular momentum is preferred. Then the 
asymptotic two-body wave function for particles interacting via short-range po- 
tentials behaves as 

u{r) = sm[Kr + x(k)] , (3.2) 

where x is the phase shift and hn is the relative momentum. The phase shift 
depends on the relative energy ti^K? /m and is at low energy given by the ex- 
pansion 

K cot [x{k)\ =-— + l-K^Res + 0{k^) , (3.3) 

Us ^ 

where a« is the s-wave scattering length and i?efF is the effective range. The 
convention applied here is that for a purely repulsive interaction the scattering 
length is positive, while for a purely attractive interaction without any bound 
states the scattering length is negative. The effective range and higher-order 
terms can be neglected at sufficiently low relative energy. Thus, at low energy 
the properties of the two-body system are basically determined by the scattering 
length ttg. 
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The s-wave scattering length for a given two-body potential V{r) can be 
obtained by solving the radial Schrodinger equation for two identical particles 
of mass m for zero angular momentum, zero energy {k = 0), and boundary 
condition u{0) = 0: 



w(r) = . (3.4) 



Outside the two-body potential the solution is a straight hne. According to a 
Taylor expansion of equation (3.2), the wave function for smah k is 

u{r) ~ {l + Kr cot [x('«)] } sin [x(k)] oc 1 . (3.5) 

Thus, the scattering length can be determined by the intersection of the asymp- 
totic wave function with zero, that is u{as) = 0. 

It is often convenient to also define the parameter Ob by 

which is the Born approximation to the scattering length Og. The last equality 
holds for a central potential. The strength of the interaction is then proportional 
to Gb- 

Since the finer details of the interaction potential are superfiuous, a finite- 
range Gaussian potential 

= ''« = ^. <") 

see dotted hne in figure 3.1, is sufficient for a study of the dependence on the 
scattering length and possibly a few more of the low-energy parameters in the 
expansion of the phase shift. The strength Vq is then related to Ob as indicated. 

Figure 3.2a shows function of the strength parameter Ob for the 

Gaussian potential. When the parameter 0,3 decreases from zero to negative 
values, the scattering length varies slowly and roughly linearly with ob for small 
Ob , until a value ^ where as diverges as a signal of the appearance of the first 
two-body bound state. For increasing attraction Us turns positive when this 
state is slightly bound. Then the scattering length decreases and turns negative 
again. This pattern repeats itself as the second bound state appears, and so on 
at each subsequent threshold. 

For a square- well potential Vsv,{r) = Vs^y,oQ{r < h) the threshold value of 

differs from the value for the Gaussian potential, but a^/aB as a function of 
ob/ob^ results in virtuahy the same curves, see figure 3.2b. This indicates that 
for simple potentials the behaviour is approximately independent of the shape. 

Table 3.1 shows the scattering length for different potential strengths Ob 
for the Gaussian potential, primarily for the cases studied in this work where 
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Figure 3.2: a) Scattoring length divided by the potential range 6 as a function 
of flB divided by b for the Gaussian potential from equation (3.7). b) Scattering 
length as divided by a-^ as a function of divided by defined as the 

value of flB where the first bound state occurs. Results are shown for the 
Gaussian potential with ag^/& = —1.1893 and for the square- well potential 
V,„{r) = Kw,oe(r < b) with a^°^ /b = -0.8225. 



I Ob 1/6 is close to unity. The Born approximation equals the correct scattering 
length only in the limit of weak attraction where the magnitude of the scattering 
length Us is much smaller than the interaction range b. 

To exemplify, in experimental work *^Rb atoms have a scattering length of 
tts — 100 a.u.* Assuming an interaction range around b = 1 nm we obtain 
tts/b = 5.29. This can be modelled by a Gaussian two-body interaction with 
ttB/b ~ —1.5, where the lowest solution corresponds to two-body bound states 
and the next accounts for the properties of the dilute gas. However, by applying 
an external magnetic field it is possible to change the internal energy levels in 
alkali atoms, e.g. in ®^Rb [CCR+00], and thereby change the scattering length 
to almost any desired value. This allows experimental studies of a large range 
of scattering lengths. 

*The scattering lengths for relevant spin states of *^Rb atoms are according to Pethick 
and Smith [PS02] all close to 100 a.u., where 1 a.u. = 0.529 • lO-^" m. 
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+1.00 
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+1.00 
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-0.3560 
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-1.1893 
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-0.84 







-1.2028 


+100 
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-0.551 


-1.00 







-1.220 


+44.5 
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-1.00 


-5.98 







-1.3380 


+ 10.0 


1 


-1.069 


-10.0 







-1.35 


+9.32 


1 


-1.110 


-15.7 







-1.50 


+5.31 


1 


-1.1761 


-100 







-6.868 


-1.00 


1 


-1.1860 


-401 







-7.6612 


-10.0 


1 



Table 3.1: The scattoring length as in units of b for various strengths of a 
Gaussian potential measured as UB/b. The number Ab is the number of bound 
two-body states. 

The short-range two-body interaction with s-wave scattering length Os has in 
mean-field contexts [DGPS99], i.e. with a Hartree ansatz as in equation (2.31), 
been modelled by the three-dimensional zero-range potential 

Vs{r) = ^^6ir) , (3.8) 
m 

where S is the Dirac delta function. Only the scattering length enters as the 
parameter characterizing the two-body interaction. This is usuahy assumed to 
be succcsful when njosp <C 1, where n is the density of the system. For this 
zero-range interaction equation (3.6) yields ob = fls, which is rarely the case for 
finite-range interactions, as is obvious for the cases illustrated in figure 3.2. 

The finite-range Gaussian potential from equation (3.7) is used in the fol- 
lowing calculations. In order to test the dependence on the short-range details 
of the interaction, a linear combination of different Gaussians was also used in 
some cases, although these results are not shown here. 

3.2 Analytical angular properties 

Before solving numerically we investigate various limits analytically. In the 
non- or weakly- interacting limit, the kinetic-energy eigenfunctions are relevant 
for understanding the properties of the many-body system. When a two-body 
bound state is present, there is a signature of it in the angular spectrum, which 
can also be studied analytically. A zero-range treatment incorporates the well- 
known asymptotic two-body behaviour into the many-body wave function. This 
leads to an equation which has an analytic solution for a very dilute system. 
Finally, we average the interactions in a way that resembles the mean field, 
i.e. all correlations are neglected. These different analytic approaches provide 
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a basis for understanding the numerical solutions, which we turn to in section 
3.3. 



3.2.1 Kinetic-energy eigenfunctions 

First, non-interacting particles, that is z; = 0, are considered. With the trans- 
formation in equation (2.20), equation (2.75) becomes 



— + ^ 9 tan^ a A 

da^ 4 2 



4>{a) = . (3.9) 



Here (^(a) is a reduced angular wave function 

4>{a) = sinacos^^^-^)/^ a </.(«) , (3.10) 

in analogy to the transformation from radial to reduced radial wave function 
for the two-body problem. Since (j) for a physical state cannot diverge at a = 
or a = 7r/2, the boundary condition for the reduced angular wave function is 
^(0) = = 0. 

Non-reduced solutions to equation (3.9) are given by the Jacobi polynomials 
V [AS65] as 

0^(a) = p[V2.(3^-8)/2l(cos2a) . (3.11) 

See further details in appendix C. The hyperspherical quantum number K is 
given by K = 2v = 0,2,4,... and denotes the angular kinetic-energy eigen- 
function with v nodes in a space. The corresponding angular eigenvalues are 

= K{K + 3N — 5). This notation is consistent with the general hypersphe- 
rical harmonics from equations (2.48) and (2.50). The lowest eigenvalue is zero 
corresponding to a constant eigenfunction Vq = I. 

Figure 3.3a shows the reduced angular kinetic-energy eigenfunctions for A'' = 
100 and the lowest three eigenvalues. The constant wave function (f>K=Q is in 
the figure represented by (t>aioi) = sinacos'^^"'^^/^ a, where |0oP then is the 
volume element in a space. The oscillations are located at relatively small a 
values. As seen in figure 3.3b, the location of the maximum changes as 1/Vn 
due to the centrifugal barrier proportional to tan^ a in equation (3.9). Thus, 
as N increases, the probability becomes increasingly concentrated in a smaller 
and smaller region of a space around a = 0. 

Some solutions may be spurious, i.e. each component </> is non- vanishing, but 
the full wave function $ in equation (2.64) is identically zero: 

N 



Equation (2.69) shows that such a component 4> with zero sum is an eigenfunc- 
tion of the angular kinetic-energy operator. Here the K = 2 eigenfunction from 
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Figure 3.3: The reduced angular wave function (pK, defined in equations (3.10) 
and (3.11), for a) N = 100 and = 0, 2, 4 and b) X = and iV = 10, 100, 1000. 
The normahzation is J^^^ da \4>K{ct)\'^ = 1- 



equation (3.11) has a vanishing angular average, i.e. 

N 

dr J2'l>K=2{aij) = , (3.13) 

see appendix C. This criterion is not identical to equation (3.12), but functions 
(j) that obey equation (3.13) are nevertheless inert to the interaction potential 
as it occurs in the Faddeev-Hke equation (2.72) and in the angular variational 
equation (2.84). Solutions Hke (j)K=2 obtained by solving equations (2.72) or 
(2.84) are therefore independent of the interactions and the eigenvalue is inde- 
pendent of p. Since the K = 2 function is spurious in this sense, it must be 
avoided when obtaining the solutions. 



3.2.2 Asymptotic spectrum for two-body states 

For large values of p, the short-range two-body potential v with range b is non- 
vanishing only when a is smaller than a few times b/p. For larger values of a, 
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the "potential-rotation" terms vRcp in the angular Faddeev-like equation (2.75) 
can therefore be omitted. 



We first assume that the rotation terms R<?f) can be neglected for smaller a. 
For Q ^ 1 substitution of r ^ \/2pa instead of a in equation (2.75) then leads 
to the two-body equation with energy 




+ V{r) - E^^^ u{r) = , 



(3.14) 



where 2mp'^E^'^'^ /h^ = X+9N/2-9 and u(\/2psinQ!) = ^(a). A two-body bound 
state with E^'^^ < corresponds to an eigenvalue A diverging towards — oo as 
— p^. Moreover, the wave function will be concentrated around r ~ 6, which 
in terms of a means a ^ 6//0 <C 1. Such solutions do not produce significant 
rotation terms, which is consistent with the omission in the derivation. The 
structure of the A^-body system is given by the fully symmetrized wave function 
for two particles in the bound state, while all other particles are far away, thus 
producing the large average distance. 

A solution to equation (2.75) that does not correspond to a two-body bound 
state has a wave function distributed over larger regions of a space. As the 
potentials then vanish for large p, we are left with the free solutions, i.e. the free 
spectrum of non-negative A values is obtained in this Hmit of large p. 

A two-body state with energy slightly below zero forces A to diverge slowly 
as — p^. On the other hand, if the two-body system is shghtly unbound, A 
instead converges slowly to zero which is the lowest eigenvalue of the free solu- 
tions. Precisely at the threshold, it seems that A should not be able to decide 
and therefore must remain constant. Thus, for infinitely large two-body s-wave 
scattering length we are led to expect that one angular eigenvalue approaches 
a negative constant for large p. Similar predicted behaviours have been con- 
firmed for three particles [NFJGOl]. In section 3.3.2 we turn to the numerical 
verification for > 3. 

3.2.3 Zero-range approximation 

A zero-range treatment of three-body systems leads to an equation which can be 
easily solved for the angular eigenvalue [NFJGOl, FJOlb]. The basic assumption 
is that when the liyperradius is large compared to the range of the interaction, 
two interacting particles in the many-body system consider each other as point 
particles. Therefore, the details of the interaction potential can be replaced by 
a boundary condition at zero separation. Moreover, it is seen from the Faddeev- 
like equation (2.75) that outside the potential range the angular equation is just 
the kinetic-energy eigenvalue equation with the solutions from section 3.2.1. 
Combination of these observations leads to an analytic solution as follows. 

The many-body wave function at small two-particle separation approaches 
the two-body wave function. The two-body wave function at low energy accord- 
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ing to equation (3.5) then behaves as 

1 du{r) 



u{r) dr 



r=0 



1 



(3.15) 



In order to compare consistently, we use the many-body wave function including 
the volume element in the angle a where a is related to the two-body distance 
r by r = \/2psma, i.e. wc assume that at small separations the two-body wave 
function u{r) is represented by a^{p,fl). Equation (3.15) then becomes 



d[a'^{p,n)] 



da 



V2p 



a=0 



(3.16) 



a=0 



$(p, n) = ^(a) + 2{N - 2)R[^-'U{a) + ^{N - 2){N - 3)R(r'V(a) ■ (3.17) 



By averaging over all angular coordinates except a we obtain 

1 

Outside the diagonal potential v{a) the solutions and eigenvalues, with 
proper boundary condition at a = 7r/2, are 

0,(a) = cos2a) , P^x) = Pp^-8)/2.V2](^) ^ (g ^g) 

X = 2u{2u + 3N -5) . (3.19) 

Since we do not restrict <pv for a — > 0, non-integer values of u are allowed. For 
small a the solutions behave as [NFJGOl], see also appendix C, 

, . ^ ^ „ . sin(^;y)r(i/+^i^) 



B = cos(7rJ^) 



A = 

2 r(z. + |) 



0Fr(:/ + i) ■ 

Then at the edge of the zero-range potential we get 



a<f{p,n) 

a=0 

d[a^p, n) 



A 



a=0 



da 



= B + 2{N-2)r[^-^^M0) 



a=0 



+ l(iV_2)(Ar-3)R(r'V.(0), 



R34 \ 



..(0) 



2r(^)r(. + |) .^p^ 



3N-6^ 



^13 



0.(0) 

1/2 



2r 



^/H^(^)V3y' 



2\ 



(3]V-8)/2 



(37V-ll)/2 ^ 



1 for TV > 3 , 



(i/+l)V3 



(3.20) 
(3.21) 

(3.22) 

(3.23) 
(3.24) 

(3.25) 
(3.26) 
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Combination of these results leads to p/ag as a function of v. 

p _ x/2r(. + |)r(3^)r(.+ 3^) ^ 

as sin(7r!y) y{u+^^)^ 

cos(7r^)r(^+^) .. Rir^V.(O) , (iV-2)(7V-3) 

At small <C 1, the square bracket yields N{N — l)/2, and then becomes 

2V2^ r(3^)7- 

The angular eigenvalue A from equation (3.19) is then 

-jV(jV-l)^A_^_i. (3.29) 

This derivation is valid when <C 1, or equivalently when p N^^^^\as\. 

As we shall see in the following section, the result in equation (3.29) can be 
obtained otherwise. However, when the treatment of equation (3.27) is numer- 
ically extended to smaller hyperradii, unmistakably wrong results are encoun- 
tered. Whether this is reminiscent of the initially expected deficiencies of the 
Faddocv-likc equation or it is a mistake in the treatment of equation (3.27) is 
presently not sorted out. 



3.2.4 Average, non-correlated effects of interactions 

As discussed in section 2.3.1, a mean-field wave function corresponds to a con- 
stant angular wave function where no correlations are included. With a non- 
correlated, constant angular wave function ^k=o = '^i'<:j 4'K=o{aij), the ex- 
pectation value of the angular Hamiltonian Hq becomes 

N 

Aif=o = {^K=o\hQ\^K=o)n = (^K=o\^Vki\^K=o) , (3.30) 

k<l 

without contribution from angular kinetic energy. Proceeding in the manner of 
the mean field we then have to assume the same ansatz for the two-body inter- 
action, i.e. the S function from equation (3.8). With this zero-range interaction, 
equation (3.30) becomes 

Bohn et al. [BEG98] did a similar calculation, but since they did not separate 
out the centre-of-mass motion, the present result for is effectively that of 
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[BEG98] with TV replaced by A/' - 1 in the T function. For TV > 1 the results 
are identical. 

We note that the angular potential from equation (3.31) coincides with equa- 
tion (3.29), i.e. the large-hyperradii derivation from the zero-range model in sec- 
tion 3.2.3. This indicates that the structure of the two-body correlated ansatz 
for the many-body wave function catches the essential information in agree- 
ment with the low-density result, equation (3.31), which corresponds to the 
mean field. 

Thus, the zero-range interaction from equation (3.8) leads to reasonable 
energies in the dilute limit. However, at larger densities (smaller p) a nega- 
tive scattering length potentially leads to unphysical behaviours. We can 
understand this problem by putting a tts/p into the radial potential, equa- 
tion (2.29), which yields a term risl that diverges faster than other terms as 
p ^ 0. We return to this problem in chapter 6. 

Thus, a zero-range two-body interaction in mean-field computations can 
lead to a collapse. This problem is not present for finite-range interactions, 
and the present method allows the use of strongly attractive potentials. The 6 
interaction furthermore does not allow a study of short-range properties such 
as bound two-body systems and similar clusterizations. Both problems arc 
overcome by using a finite-range potential in the present model. When p is 
much larger than the potential range h, the expectation value of a finite-range 
potential is of the same form as A^ in equation (3.31) 



with the Born approximation ae from equation (3.6) instead of the real scatte- 
ring length ttg. 

In the opposite limit, when p <^ b, the result is strongly dependent on the 
shape of the potential. For example, the Gaussian potential from equation (3.7) 
yields 



As seen from these two limits there are some scaling properties for finite-range 
potentials. The angular eigenvalue at a given N value depends only on ae/ft 
and p/b. For a Gaussian potential we have 




(3.32) 




(3.33) 



Vkl = 



2mp^Vo ^ Sob /P\ 2 2(^/6)2 si„2„^ 



^b\bJ 



(3.34) 



which implies that for a given value of Ob/6, the angular eigenvalue A is only a 
function of p/b. The radial potential U from equation (2.29), which we return 
to in chapter 4, can be scaled as 
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where b^ = ^h/{mijj) is the characteristic length for a harmonic trap of an- 
gular frequency lo. The scaled energy 2mb'^E/h^ is then for a given N value 
only a function of ae/?) and bt/b. These scaling properties are useful in model 
calculations. 

3.3 Numerical angular solutions 

In the previous section we discussed solutions to the angular equation in the 
presence of no interactions, in the case of two-body bound states, and in the 
zero-range limit. However, solutions with general two-body interactions have 
to be obtained numerically, which is the quest of the present section. We first 
comment on the numerical procedure before discussing properties of the angular 
eigenvalues and wave functions. 

3.3.1 Numerical method 

The angular eigenvalue equation was rewritten in chapter 2 by a variational 
technique as the second-order integro-differential equation (2.84) in the variable 
a, where ri2 = v^psina. For neutral atoms in recent trapping experiments 
the interaction range is very short compared to the spatial extension of the N- 
body system. Then this equation simplifies to contain at most one-dimensional 
integrals. The validity of the approximations only relies on the small range b of 
the potential, whereas the scattering length as can be as large as desired. 

Even though the complexity of the angular equation does not increase as the 
number of particles increases, the numerical solutions become harder to handle 
for large N. The origin of this problem is the sharp peak in the angular volume 
element for large N, see section 3.2.1. 

Expansion on kinetic-energy eigenfunctions 

A usual method within the hyperspherical formalism is to expand the angular 
wave function on kinetic-energy eigenfunctions [Lin95, BEG98]. Such an ex- 
pansion is successful when the physical extension of the system is comparable 
to the interaction range. The hyperspherical harmonics contain oscillations at 
angles of the order of magnitude a ~ 0{1/K), so for a given hyperradius we 
need K values of the order of Kma.x ^ ^{p/^) to describe potentials limited 
to a < b/ p. Thus, the angular kinetic-energy eigenfunctions constitute an in- 
effective basis at large hyperradii since the diagonal potential in this case will 
be sharply peaked around q = 0, and a huge number of terms is necessary to 
account for the correct behaviour of the wave function around a = 0. 

For trapped particles the scale of the system is determined by the trap length 
6t which for atomic gases usually is of order /im. Since the interaction range 
b usually is in the nm region, an expansion on kinetic-energy eigenfunctions 
converges slowly and is not appropriate for the present treatment. 
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Finite differences 

Instead of an expansion on hypcrsphcrical harmonics we choose a basis of dis- 
crete mesh points distributed in a space ^(a) (p = [(j){ai), . . . ,4){aM)] to 
take into account the short range of the potential and to keep sufficient infor- 
mation about small a. Derivatives are then written as finite differences [KM90] 
and integrations like R(/)(a) of equation (2.88) can be expressed in matrix form, 
i.e. R^(q;) R(f). 

Numerical computation of the integrals becomes increasingly difficult with 
decreasing interaction range. This can be understood in terms of the a coordi- 
nate, since the potential at a given p and a given range b of the interaction, is 
confined to an a region of size Aa ~ b/p, which for Bose-Einstein condensates 
easily becomes very small and thus cannot be handled directly numerically. 

Recently the method of finite elements was applied to the Faddeev-like equa- 
tion (2.75). With finite elements the basis functions are smooth and yield more 
reliable matrix elements, especially those involving derivatives due to the kinetic 
energy. This proves easier to handle, but is presently not implemented for the 
angular variational equation (2.84). For details about finite-elements methods 
see references in Press et al. [PFTV89]. 

Unless stated otherwise, the following numerical results are obtained with 
the method of finite differences. 

3.3.2 Angular potentials 

The angular eigenvalue depends on the number of particles, on the size of the sys- 
tem through the hyperradius, and on the two-body potential. Figure 3.4 shows 
the angular eigenvalue for the particle number N = 20 and various Gaussian 
potential strengths. Only the lowest Aq is shown unless otherwise indicated. 

The long-dashed line shows the calculation for a purely repulsive interaction 
with positive scattering length. Here the angular eigenvalue approaches zero at 
large hyperradii approximately as The thin, soHd fine shows oc Us/p from 
equation (3.31) for the same scattering length. These two curves almost coincide 
at large hyperradii. The short-dashed curve shows the angular eigenvalue for a 
slightly attractive two-body interaction without any two-body bound state and 
with negative scattering length. This angular potential approaches zero from 
below as 1/p, also in agreement with equation (3.31). For a larger attraction, 
when the scattering length becomes very large, the angular eigenvalue (thick, 
soHd line) is almost constant for a large region of hyperradii. This agrees with 
the expectations in section 3.2.2. For a slightly larger attraction the scattering 
length turns positive and a two-body bound state forms. Then (dot-dashed line) 
the lowest angular eigenvalue at some point diverges to minus infinity. For even 
larger attraction the binding energy of the bound state increases and A diverges 
faster, see the sequence of the dot-dashed, double-dashed, and triple-dashed 
lines. 

The dotted line shows the angular eigenvalue for the next angular solution 
for the strongest attraction. This approaches zero from above as 1/p, which 
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Figure 3.4: Angular eigenvalues for N = 20 and parameters {aB/b,as/b) as 
shown on the figure. A star refers to the first excited state. For ae/fe = —1.1893 
we have Ug/b = —85601, see table 3.1, which here is denoted by — oo. 



resembles the behaviour for a purely repulsive interaction (long-dashed line). 
This illustrates the use of the terms "effectively repulsive" or "effectively attrac- 
tive" in the mean field, depending on the sign of as even though the interaction 
potential might be purely attractive. See related comments by Geltman and 
Bambini [GBOl]. 

For one hundred particles figure 3.5 shows the lowest angular potential for 
various attractive interactions. Qualitatively the same behaviours as for = 20 
are observed. When Us = —b (sohd line) the system has no bound two-body 
states. The lowest angular eigenvalue is zero at p = 0, decreases then through 
a minimum as a function of p, and approaches zero at large hypcrradii as as/ p. 
A larger attraction (broken lines) decreases all angular eigenvalues for all p 
values. The details at smaller hyperradii hardly change with large variations of 
the scattering length. However, at larger distances the approach towards zero 
is converted into a parabolic divergence as soon as the scattering length jumps 
from negative (dotted line) to positive (dot-dashed line) corresponding to the 
appearance of a bound two-body state. The faster divergence (double-dashed 
line) is again observed for increasing binding energy. 

The characteristic feature for both cases N = 20, 100 is the large-distance 
asymptotic behaviours. For repulsive potentials all eigenvalues are positive and 
the lowest approaches zero from above. The higher eigenvalues would then 
converge to K{K + 3N-5) as 1/p, where K = 4,6,8.... The solution for K = 2 
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Figure 3.5: The lowest angular eigenvalues A for iV = 100 bosons interacting via 
a Gaussian two-body potential V{r) — Vb exp(— with zero or one bound 
two-body states. The scattering lengths Us/b are indicated on the figure. 



is not allowed, corresponding to removal of the non-physical spurious solution, 
see section 3.2.1. 

For weak attractions the lowest A is negative and approaches zero from be- 
low as 1/p. The higher angular eigenvalues approach, again, K{K + 3N — 5) 
corresponding to the spectrum for free particles. The constant of proportional- 
ity to p^^ for the lowest eigenvalue is qualitatively recovered as the predicted 
dependence on Ug. Calculations with a two-body potential as a linear combina- 
tion of different Gaussians (not shown) confirm that the large-distance angular 
potential only depends on the scattering length as as in Xs. 

In the presence of a two-body bound state the divergence as —p^ refiects 
the corresponding two-body binding energy, see equation (3.14). Generally, an 
attractive finite-range interaction can support a certain number Ab of two-body 
bound states for both positive and negative scattering lengths. Then the lowest 
angular eigenvalues, Ao, Ai, . . . , Atv'b-Ii describe these bound two-body states 
within the many-body system at large hyperradii, i.e. they diverge to — oo as 
seen in figure 3.5. 

The next eigenvalue Atv^ converges to zero at large distance and corresponds 
to the first "two-body-unbound" mode. The higher eigenvalues would then, once 
more, converge to K{K + 3N — 5). Increasing the attraction to allow another 
bound two-body state would then shift the asymptotic spectrum such that one 
more eigenvalue diverges while the non-negative energy spectrum remains un- 
changed. This yields qualitatively the same asymptotic spectrum for the un- 
bound modes irrespective of the number of bound states below. This invaluably 
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eases the computations, i.e. all the bound states of the two-body system are 
not needed in order to describe the unbound modes of the many-body system. 
Therefore, the two-body interaction does not have to be the real two-body inter- 
action, which allows all the two-body bound states that are known to exist, but 
the interaction potential can be written in a way that accounts for the investi- 
gated properties. This is the case for the potential from Geltman et al. [GBOl], 
equation (3.1), and also for the Gaussian potential, equation (3.7), applied in 
the present work. 

These properties of the two lowest eigenvalues in the presence of one two- 
body bound state are evident in figure 3.6. The lowest eigenvalue (dashed 
curve) diverges to minus infinity proportional to . This corresponds to the 
bound state. The second eigenvalue (solid curve) is negative at smah hyperradn, 
but turns positive at larger and approaches the asymptotic behaviour of \s oc 
as/p (dotted curve, see details in the inset). Since the second eigenvalue at 




Figure 3.6: The two lowest angular eigenvalues (dashed and solid curves) for 
N = 100, cis/b = +10, and one bound two-body state. The dotted curve is Xs 
for the same scattering length. 



small and intermediate hypcrradii is negative, this might allow a self-bound 
system located at distances far inside and independent of a confining external 
trap potential. This feature is absent in a description with overah repulsive 
potentials, corresponding to positive scattering lengths, for example the zero- 
range interaction with > 0. Then no attractive part is possible. 

At each threshold for the appearance of a new bound two-body state, one 
eigenvalue asymptotically approaches a negative constant as in figure 3.4. This 
eigenvalue is responsible for the structure of the A^-body system for very large 
scattering lengths. This refiects the transition from unbound to bound two- 
body states, that is the transition from convergence towards zero as — 1/p to 
divergence as —p^, see section 3.2.2. 

We finish the discussion of the angular eigenvalue from the variational equa- 
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tion (2.84) by comparing with the results from a finite-element treatment of the 
Faddeev-like equation (2.75). Figure 3.7 shows the results from the Faddeev-like 
equation for the same parameters as in figure 3.4. At large hyperradius the re- 
sults agree, whereas they differ as p ^ 0. This is probably due to the short-range 
approximation of section 2.4.3, although we recall the non-variational nature of 
the Faddeev-Hke equation, as discussed in section 2.4.1, as another possible 
source. However, due to the importance of correlations higher than two-body in 
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Figure 3.7: Angular eigenvalues for = 20 and parameters {aB/b,as/b) as 
shown on the figure obtained from the Faddeev-like equation (2.75). 



the denser regions, corrections to these short-distance results would be in order 
even without the short-range approximation. Moreover, for a description of a 
dilute many-boson system we do not need the details at such short distances, 
so they are not considered in the following. 



3.3.3 Angular wave function 

The total angular wave function is determined as the sum of two-body compo- 
nents in equation (2.64). Figure 3.8 shows the lowest component wave function, 
reduced as in equation (3.10), for a two-body potential with one bound two-body 
state. With increasing p the amplitude concentrates at smaller and smaller val- 
ues of a. This refiects the convergence towards the two-body bound state in 
agreement with the transformation ri2 = \/2psiuo, see section 3.2.2. The nu- 
merical recovery of this behaviour is essential, since otherwise the large-distance 
properties cannot be described. 
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Figure 3.8: The lowest reduced angular wave functions for iV = 20 and = 
— 1.506, as = 5.315 for three values of the hyperradius. This potential has one 
bound two-body state. 



The angular eigenfunction varies with the strength of the interaction. Ex- 
amples of this variation are shown in figure 3.9a. The lowest non-interacting 
wave function (thin, solid line) has only nodes at the endpoints. The repulsive 
case shows an oscillation (dashed line) which lowers the angular energy due to 
the rotation terms. The fast change at small a, which is emphasized in fig- 
ure 3.9b, is typical for interacting particles. The wave function for the excited 
state (dotted line) has an additional node. The corresponding lower-lying wave 
function was shown as the dashed line in figure 3.8. 

The wave function for infinite scattering length (thick, soHd fine in figure 3.9) 
corresponds to an interaction where the two-body bound state is at the threshold 
for occurrence. This eigenfunction resembles those where a bound two-body 
state is present, compare with the results shown in figure 3.8. However, now 
(thick curve of figure 3.9b) the wave function is located at larger a values. 

The properties of the component of the angular wave function is further 
illustrated by the second moment defined by 

{rl2)4> = 2p2(0i2| sin^ a|(/.i2) . (3.36) 

A number of these moments for different interactions are shown in figure 3.10 
as functions of p. For states obtained from repulsive potentials, moderately 
attractive potentials without bound two-body states, and for excited states of 
positive A, the moment (rij)^ increases proportional to for large p. This 
resembles the behaviour of the expectation value in the lowest angular state 
for a non-interacting system, i.e. K = 0, where (^"12)0 = 2/3^/(A^ — 1). The 
qualitative explanation is that large p implies the limit of a non-interacting 
spectrum with the corresponding non-correlated wave functions. 

In contrast, a difi^erent behaviour is observed when the potential can bind 
two particles, i.e. (^^2)0 approaches a constant at large p. The angular equation 
in this limit approaches the two-body equation (3.14). The wave function in 
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Figure 3.9: a) Angular wave functions for N = 20 and p = 5006 for different 

interaction parameters {a-Q/b,a.s/h) as shown on the figure. The K = {) curve 
corresponds to a non-interacting system. A star refers to the first excited state, 
b) The same as a), but with logarithmic a axis. 



the zero-range limit converges to u{r) — exp(— r/os). The second moment is 
then found as {u\r'^\u) = a1/2, which in the limit of large p reproduce the 
constant values for (rfj)?; when a.,/6 = 9.32 and ag/h = 5.31, i.e. the double- 
and triple-dashed hues in figure 3.10 approach 9.32^/2 ~ 43 and 5.31^/2 ~ 14, 
respectively. 

Expressed differently, when a two-body bound state is present, the angular 
wave function is at increasing p squeezed inside the potential since the range 
in a space decreases proportional to p~^ . This implies {4>i2\ svc? a\cf)i2) oc 1/ p^. 
The distance between a pair of particles is therefore independent of p at large 
values of p. This means that pairwise the two-body bound state is approached 
while all other particles are far away. The symmetrization does not affect this 
conclusion. Thus, apart from this symmetrization of the many-boson wave 
function, the attributes of the many-body system in the presence of this two- 
body bound state show only small deviations from the well-known properties of 
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Figure 3.10: The second moment {r\2)<p as a function of hyperradius for N = 
20 for solutions to tiie angular variational equation with different interaction 
parameters specified in the figure by {aB/b,as/b). Also shown is the K = 
value. A star refers to the first excited state. 



the isolated two-body bound state. 

At the threshold for two-body binding, that is for infinite scattering length, 
the intermediate behaviour once again emphasizes the transition from bound to 
unbound, see the thick, solid line in figure 3.10. 

3.4 Summary 

Further numerical analysis allows us to construct a parametrization for the be- 
haviour of the lowest angular eigenvalue for attractive two-body interactions in 
two different regimes: i) no bound two-body states and Ug < 0, and ii) > 
and one bound two-body state of energy E^'^h These details are previously 
published [SF.J03a] and collected in appendix D. Here we summarize the re- 
sults, illustrate them, and then comment on them in relation to the previous 
observations. 

3.4.1 Parametrization 

For small hyperradii p < po = 0.87iV^/^(6/|as|)^/^6 we use for all the per- 
turbation result obtained as the expectation value of the two-body interaction 
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V{r) in a constant angular wave function, i.e. for N ^ 1 

mV{0)N^p'^ 



for p < po . 



(3.37) 



For hyperradii exceeding the lower limit pa the analytic expressions from equa- 
tions (D.3), (D.4), (D.9), and (D.IO) are expressed as 



\,{N,p) = -\XsiN,p)\{l + 
1 



cxp 



|Aoc(W)| 
|A^(W,p)| 
|A'^'(P)I 



P ) 
when Os < 
when a. > 



(3.38) 



for p> Pq 



with from equation (3.31) and 

Aoo(A^) = -1.59iV^/3 



m\as\ 



(3.39) 
(3.40) 



The number c approaches unity when the scattering length becomes very large. 
The factor (1 + Q.92N''/%/ p) reflects dependence on the finite range h of the 
Gaussian two-body interaction. At p ~ A''^/^|as| we find A^ ^ \^ ^ \^'^\ 

The results of the parametrizations in equations (3.37) and (3.38) are illus- 
trated in figure 3.11 for N = 100 and various scattering lengths. The pronounced 




Figure 3.11: The angular eigenvalue A, equations (3.38) and (3.37), for N = 100 
as function of p for the different scattering lengths given on the figure in units 
of the range as/b. 

deep minimum at p ^ po is in the region depending on the two-body potential 
and reflects the qualitative behaviour of the lowest angular eigenvalue. After 
this strongly attractive region at small p the eigenvalues approach zero. As the 
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size of the scattering length increases, the eigenvalue develops a plateau at a 
constant value Aqo independent of . Eventually at large p the eigenvalues van- 
ish as Xs when < and diverge to — oo when > 0. This is comparable with 
the sequence of the short-dashed, thick-sohd, and dot-dashed hnes in figure 3.4. 

When Ug < 0, the analytic and the correct eigenvalues both exceed the 
asymptotic zero-range result, i.e. Aa > Xg for all hyperradii. This means that 
the ground-state energy is higher than the energy obtained with the zero-range 
interaction. Thus, the ground state energy from the present model is higher than 
the mean-field energy. The origin of this sequence of energies is that the zero- 
range interaction inevitably leads to diverging energies for smaller distances. 
The present model avoids this non-physical short-range collapse. 

When ttg > the interaction is effectively repulsive at large hyperradii and 
an analytical expression in this case for the second angular eigenvalue obeys 
Aa < A^ for all hyperradii, due to the divergence of Xs +oo as p ^ 0. 
Correspondingly, the energies are smaller than the zero-range mean-field result 
in the positive-ttg case. 



3.4.2 At the threshold 

At intermediate hyperradii, that is when 

6<^<|a,|, (3.41) 

the angular eigenvalue as obtained from equation (3.38) is independent of both 
the short-range details of the two-body interaction and the scattering length. 
Then A approaches a constant value given by equation (3.39) as Aoo — — 1.59A^''/^. 
This plateau value can be estimated by considering the angular eigenvalue for 
a two-body bound state: 

X^^\p)c.-^. (3.42) 

The plateau terminates at a hyperradius Pa where this two-body angular po- 
tential intersects with Xs from equation (3.31), i.e. 

X^iN) = X^'Hpa) = XsiN,pa) ^.&V2^ . (3.43) 

2 V TT Pa 

Combination of equations (3.42) and (3.43) yields 

Pa ^N^'lasl, (3.44) 

X^{N) ~ -^Af^/3~-1.65iV^/3, (3.45) 

which is in agreement with the numerical results in equations (3.39) and (3.41). 
However, the N dependence cannot be predicted from the angular equation 
since the result is an interplay between the various terms in equation (2.84). 
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The symbol Aoo is chosen for this constant since the relevant p region extends 
to infinity in the limit of infinitely large scattering length. With no bound two- 
body states {as < 0) the lowest angular eigenvalue approaches zero at larger 
hyperradii, whereas it diverges towards — oo as when a bound two-body state 
is present (Ug > 0). On the threshold for a two-body bound state Ug = ±oo 
and the angular eigenvalue therefore remains constant. For finite, but large, Qg 
the eigenvalue lingers and cannot decide which way to go until the hyperradius 
exceeds a size Pa proportional to the scattering length given by equation (3.44). 

3.4.3 Discussion 

The two-body correlations built into the many-body wave function are evident in 
the properties of the angular wave function. The particles feel pairwise repelled 
or attracted to each other, which is refiected in the average two-body distance in 
the two-body amplitude. The presence of a two-body bound state is described 
by an angular adiabatic potential proportional to the two-body binding energy. 
The angular wave function in this limit equals the wave function for the two- 
body bound state. 

The angular adiabatic potential reflects the effective interaction between the 
bosons. We recovered numerically the scattering length as the determining pa- 
rameter for a dilute system with large average separation. Deviations at larger 
densities resulted in a parameter-free effective interaction Aoo which is inter- 
preted in simple physical terms as the transition between the shape-dependent 
and the scattering-length-dependent regions. The properties of the lowest an- 
gular eigenvalues are collected in table 3.2. 
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Table 3.2: The behaviour of the lowest angular eigenvalues at large hyperradii as 
a function of the number Af-g of bound two-body states and for different regions 

(2) 

of the scattering length. The attraction increases through the sequence. £"„ is 
the energy of the n'th two-body state. 

Calculations with a correlated Jastrow wave function with the right be- 
haviour at small interparticle distances and with realistic interactions [CHM+02] 
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and with various finite-range potentials [BGOl] confirm that the ground-state 
energy of a dilute boson system only depends on the scattering length and not 
on the details of the potential. The advantage of the present model is that it re- 
sults in a relatively simple one-dimensional differential equation which provides 
analytical results in some Hmits, i.e. the scattering-length-only behaviour is re- 
covered analytically in section 3.2.3. Furthermore, from the two-body ansatz 
for the wave function no further assumptions arc necessary in order to obtain 
the large-distance scattering-length-only signatures. 



Chapter IV 



Hyperradial confinement and condensates 



In chapter 3 we studied the attributes of the solutions to the iiyperangular 
equation for a fixed value of the hyperradius. This frcczcd the variation in 
the average distance between the particles, but nevertheless showed a range of 
characteristics depending on the nature of the two-body interaction. In this 
chapter we complete the treatment of the degrees of freedom in the centre-of- 
mass system by studying the radial equation and the properties of its solutions. 

Besides the contributions from kinetic energy and interactions, the radial 
equation contains a term due to an external field acting on the particles. As 
shown in chapter 2 this separates for a harmonic field nicely into a centre-of- 
mass part and a hyperradial part. The inclusion of such a term is discussed 
in section 4.1. Then the hyperangular contributions due to interactions are 
included in section 4.2 in a study of the properties of the full radial potential 
and the solutions to the radial equation. Section 4.3 presents more details about 
negative-energy states, which include the Efimov-like states that are described 
further in section 4.4. In section 4.5 we discuss condensation before summing 
up in section 4.6. 

4.1 Trapped bosons 

In experiments neutral atoms, for instance evaporated sodium atoms [DMA+95], 
are cooled and trapped by lasers, and then held and further cooled in magnetic 
fields which interact with the magnetic moments of the particles. In a common 
set-up, the time-averaged orbiting potential, a static magnetic field is combined 
with a rotating magnetic field [PS02]. This effectively generates a harmonic- 
oscillator field in which all particles move, e.g. for particle i we have 

V^trap(T-i) = ^rn{ulx^i + ulvl + i^lzl) , (4.1) 

where the position of particle i is ri = {xi,yi,Zi), and the angular frequencies 
along the coordinate directions q = x,y,z are denoted by u)q. These angular fre- 
quencies u)q depend on the magnetic moments of the atoms and on the strengths 
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of both the stationary field and the time-varying magnetic field. Experimen- 
tally it is possible to obtain the same effective frequency along two and 
y, and a different frequency along the third axis, z. E.g. for ^^Rb atoms the 
effective frequencies ^'q = Wg/(27r) in recent experiments are Vx = = 17.5 Hz, 
and i^z ~ i^x/^ [DCC"'"01]. In terms of the trap lengths bq = ^/h/JmoJq) this is 
bx =by = 2591 nm and bz ~ V^bx- 

We will address the general geometry in chapter 7 and here restrict ourselves 
to a spherically symmetric field, oj = uJx = = ij^z, which leaves a central 
potential 

Krap(ri) = ^mwV^ . (4.2) 

Put differently, we treat the axial field as spherical with oj = ^uj^ujyUJz as 
the geometric mean angular frequency. A set of parameters which we will use 
frequently is for ^'^Rb-atoms with oscillator frequency I'trap = uj/{2tt) = 200 Hz 
[BEG98], thus yielding bt = ^/Kf{ muj) = 763 nm. All lengths are then scaled 
in units of the typical interaction range 6 ~ 10 a.u., which leads to bi/b ~ 1442. 

In the case of the free angular solutions from section 3.2.1, we have \k = 
K{K + 3A^ — 5) with X = 0, 2, 4, . . .. In the general case, i.e. when we include 
dependences beyond the s-waves in one hyperangle, we can replace K with 
Kn-1 from equation (2.49). The radial solutions are analytically obtained from 
equation (2.30) with the radial potential from equation (2.29) with zero couphng 
terms. The radial wave functions are then given by 

/„(p) = e-.V(26D^..,.+i£a.,.+i/2) i^^^ ^ 37V_-6 ^ ^ ^ ^^^^ 

where £„ for n = 0, 1, ... is the generalized Laguerre polynomial with n radial 
nodes. Here In,k plays the role of a generalized angular momentum due to 
the kinetic energy of the many-body system. Especially, Z2,o = reproduces 
the famihar behaviour of the harmonic-oscillator solutions for the two-particle 
system [BJ83]. The energy is = hw^N -l)/2 + 2n+K] with the subtraction 
of the centre-of-mass ground-state energy ihui/2. 

The ground state with n = is usually associated with the mode of Bose- 
Einstein condensation, see mean-field approaches [PS02, PS03] or related hyper- 
spherical approaches [BEG98, WM99]. In terms of temperature the condition for 
onset of Bose-Einstein condensation is that the thermal length scale given by 
k^T Ji? /{ml\) is larger than the average distance f between particles [PS02]. 
For bosons in a trap f \/ k^T/ (muj'^) /N^^^ , which yields that condensation 
occurs when 

kBT < N^l^huj . (4.4) 

When this equation is fulfilled, a large number of atoms prefer the ground 
state which is a signature of the condensate. The criterion is equivalent to 
a sufficiently large level spacing of the modes of the harmonic oscillator. In 
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the present thesis we assume no finite-temperature effects, corresponding to a 
sufficiently large level spacing. 

In this non-interacting picture all quantum modes of the many-boson sys- 
tem arc represented by single-particle levels. This means that the Bosc-Einstcin 
distribution can be applied to the non-interacting single-particle levels. In ge- 
neral the interactions complicate matters. In the following sections we discuss 
the structure of the radial solutions when interactions are included, and then in 
section 4.5 return to this problem. 

4.2 Radial potential and solutions 

The effect of interactions for a fixed value of the liyperradius was discussed in 
chapter 3 in terms of the angular potentials. These in themselves tell only a 
part of the story about the many-boson system. As indicated above for the 
non-interacting case, the radial equation is the next step in obtaining knowl- 
edge about the physical properties. The angular potentials and angular wave 
functions then enter the effective radial potential and transfer information about 
the interactions to quantities like energy and size of the system. 

For a dilute Bose gas the couphng terms of equation (2.28), which in the 
non-interacting case are identically zero, contribute at most about 1 % com- 
pared with other terms of the full radial equation (2.27). In the following all 
coupling terms are therefore omitted and the solutions to the uncoupled radial 
equation (2.30) are considered. This way only the angular potential A itself 
plays a role and additional information from the angular wave function ^(p, fl) 
is neglected. We should however bear in mind that coupling terms might play 
a role at larger densitites or scattering lengths. The radial potential then con- 
sists of three terms, where the repulsive centrifugal barrier and the confining 
external field both are positive. The interaction term can be either repulsive 
or attractive. The combination has structure depending on the strength of the 
interaction. 

The lowest potential for the non-interacting system, which was discussed in 
section 4.1, is shown as the thick, dashed fine in figure 4.1a for N = 20 particles. 
This has a global minimum at values of the liyperradius given by the trap length, 
that is at p ~ Ptrap = y^N/2ht ■ The schematic character of this non- interacting 
potential is representative also for a very weak two-body attraction and for a 
purely repulsive two-body potential. Corresponding solutions arc confined to the 
region between the infinitely large potential walls at smaU and large hyperradn. 

Also shown in figure 4.1a are the lowest two radial potentials for a Gauss- 
ian interaction with no bound two-body states and a small, negative scattering 
length. The deviations from the weakly-interacting case are substantial. For the 
lowest (thick, solid fine) a second minimum has developed due to the attraction 
between the bosons. This dominates at large densities, i.e. at small hyperradii. 
A barrier separates this global minimum from another minimum at large hy- 
perradii, see details in figure 4.1b. This second minimum almost coincides with 
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Figure 4.1: a) Radial potentials Uq (thick, solid line) and Ui (dotted line) from 
equation (2.29) corresponding to the lowest two angular potentials for iV = 20, 
tts/b = —0.5, and bt/b = 1442. The thick, dashed Hue shows the lowest non- 
interacting potential, that is with A = 0. Horizontal Hues show the lowest two 
energy levels in Uq. b) Details at larger hyperradii with the next nine energy 
levels in Uo- Here Uq and the curve for non-interacting particles are hardly 
distinguishable. 



the minimum for the non-interacting case and these are hardly distinguishable 
in the figure. 

With this potential the diagonal radial equation is solved. The solutions 
can be divided into groups related to either the first or the second minimum. 
The lowest two radial eigenstates in the lowest potential have negative energies, 
indicated as horizontal lines in figure 4.1a, and the hyperradial wave function 
is located in the global minimum at relatively low hyperradii. They are truly 
bound states as they cannot decay into continuum states at large hyperradii. 
Their properties are independent of the external trap which only has an influ- 
ence at much larger distances. These self-bound A''-body states might decay 
into lower-lying states consisting of various bound cluster states, e.g. a number 
of diatomic or triatomic clusters. We discuss this further in chapter 6. The pos- 
sibility of self-bound many-body systems even though the two- and three-body 
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subsystems are unbound is also discussed by Bulgac [Bul02], who, however, 
considers a three-body interaction strength as a determining parameter for the 
properties of the self-bound many-boson system. 

The group of states in the higher-lying minimum all have positive energies. 
These radial eigenstates are located in the trap minimum at larger hyperradii, 
see figure 4.1b, with approximately equidistant spacing as for the non-interacting 
oscillator. The lowest of these "trap states" can be interpreted as the state of 
the condensate. Thus, the structure of the "trap states" is similar for effectively 
attractive and repulsive interactions, i.e. for positive and negative scattering 
lengths. However, an attraction produces a series of lower- lying states at smaller 
hyperradii. 

The radial potential Ui corresponding to the second adiabatic potential Ai 
is shown as the dotted line in figure 4.1. This contains larger contributions from 
hyperangular kinetic energy, but still has a second minimum at small hyperradii. 
Otherwise the structure is the same with a barrier and a local minimum at larger 
hyperradii. 

Increasing A'' leaves quantitatively the same features for pure repulsion. At- 
traction leads to a decreasing barrier at intermediate hyperradius and at some 
point this barrier vanishes altogether. At the same time the attractive minimum 
at smaller hyperradius becomes deeper. This leads to an increasing number of 
bound states in this minimum as a function of N. Figure 4.2 shows Uq for 
a larger number of particles, N = 100, and doubled scattering length. The 
increased effective attraction is pronounced at large densities, that is at small 
p. The barrier height is now small compared to the potential depth at small 
hyperradii. The rather deep and narrow minimum occurs for TV = 100 about 
150 times the range of the interaction. This corresponds to a root-mean-square 
two-body distance of about 15 times the interaction range b. 

As the scattering length increases, the barrier disappears and the effective 
potential inside the trap approaches the behaviour characteristic for Efimov 
states. We return to a discussion of such states in section 4.4. 

We also study the hyperradial wave function F which tells about the root- 
mean-square displacement tr from the centre of the system defined by 

1 ^ 2 
^n^^j:ir.-Rr = '^. (4.5) 

i=l 

It is shown in figure 4.3 for various scattering lengths. The thick, soHd curve 
shows the non-interacting result. When interactions are included, the expected 
result turns up, i.e. that repulsion forces the particles away from the centre 
whereas the opposite holds for attraction. 

4.3 Self-bound many-body states 

Since the external field is negligible when p <C ^/Nbt, the radial potential is 
negative when A -|- (3A'' — 4)(3A'' — 6)/4 < and p is sufficiently small. Then 
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Figure 4.2: a) Radial potential Uq from equation (2.29) corresponding to the 
lowest angular potential for N = 100, Ug/b = —1, and bt/b — 1442. Also shown 
as horizontal lines are the negative energies Eo,n, n = 0, . . . ,46, in the lowest 
potential in the uncoupled radial equation (2.27). b) Detail at larger hyperradii. 
The energy of the first oscillator-hke state is shown as a horizontal line close to 
zero. 



self-bound many-body states with negative energies and finite extensions are 
possible. The radial equation corresponding to the relatively weak attraction 
between the twenty bosons in the potential shown in figure 4.1 has two negative- 
energy solutions with the wave function located in the global negative minimum. 

With the parametrization in equation (3.38) figure 4.4 shows the analytical 
radial potential, equation (2.29), corresponding to one of the angular eigenvalues 
from figure 3.11. The radial potential is negative in a large range of hyperradii, 
which can be divided into three different regions. For small hyperradii, region 1, 
the radial potential has a minimum. For intermediate hyperradii, Efimov region, 
the angular potential is constant and therefore the radial potential behaves as 
— 1/p^. This is from figure 3.11 seen to appear for p/b between 10^ and 10^. 
For large hyperradii, region 2, that is when p > iV^/^|as|, the angular potential 
behaves as — 1/p, so the radial potential vanishes as —1/p^. Finally the trap 
cx p^ dominates with positive contributions at large hyperradii p \fNbt. 

With the method described in [KMW02] it is possible to estimate the number 
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Figure 4.3: The probability distribution |^oP for the rms separation from the 
centre of the system fa defined by f|j = J2iLii'''i — R)'^/N = p^/N for A'' = 20, 
bt/b = 1442, and scattering lengths indicated as ttg/b. The angular potential 
was obtained from the parametrization in equation (3.38) for Og < and from 
A(5 for ttg > 0. The normalizations are different. 

N of bound states in the different regions, i.e. 

j dp^\m-){p)\, (4.6) 

where U^~^ {p) denotes the negative part of the radial potential U (p). The bound 
states in this potential can be divided into groups according to their hyperradial 
extension. The total number of such states is written as N = Ni + Ne + A/2 
where A/i, Ae, and M2 are the number of states located respectively in the 
attractive pocket at small hyperradii, in the intermediate — 1/p^ region, and at 
hyperradii large compared with the scattering length. 

The analytic expressions for the angular potential from equations (3.37) and 
(3.38) yield the crude estimate that the number of self-bound states in the 
pocket is A/i ~ X.ZN^/"^. The outer region supports bound states when the trap 
length bt is sufficiently large, that is bt ^ A/|as|, and analogously the number 
is estimated to be M2 — 0.78A/'''/^. The intermediate region is considered in the 
following section. 

4.4 Efimov-like many-body states 

When the scattering length is large, the three-body system exhibits the so- 
called Efimov effect [Efi70] where many bound three-body state turns up. In 
the following we envistigate the properties of the many-body system in this 
Efimov regime. Much of the formulation is quite similar to that in a recent 
description of three-body Efimov states [NFJGOl]. 
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Figure 4.4: Analytic radial potential obtained from equations (2.29) and (3.38) 
for N = 100, a^/b = -W^, and bt/b = 1442. 



A large scattering length implies through the eigenvalue from equation (3.38) 
an intermediate region in hypcrradius where the angular potential is almost 
constant, see section 3.4.2. More specifically, when 



b< 



NVe 



< Ids 



(4.7) 



then equation (3.39) yields A ~ Aoo = — 1.59A^''/"^ and two of the terms in the 
radial equation add to a negative value. At the threshold for binding of the 
two-body system, i.e. \as\ = oo, the radial potential in equation (2.29) then has 
the form 



U{P) 

e 



2m 



-e - 1/4 



bf 



_ iSN-A)i3N-e) _ 1 ™ ^_^g^,/3 . 
4 4 



(4.8) 
(4.9) 



This implies that no repulsive barrier is present. Then the effective potential 
behaves as — until the trap dominates. 

Figure 4.5 shows the radial potential for N = 20 and infinite scattering 
length corresponding to A^c ^ —1340 or ^ 584. Deviations from the form in 
equation (4.8) are only present at small hyperradii due to the finite range of the 
interaction. 

Without the external potential the 1/p^ potential in equation (4.8) would 
produce infinitely many radial solutions to the non-coupled radial equation (2.30). 
The radial wave function for these states would behave like 



foo{p) = v^sin 



imn(-f 

Psc 



(4.10) 



with some hyperradius scale Psc- The energies and mean-square hyperradii for 
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Figure 4.5: The lowest radial potential for A'' = 20, \as\ = oo, and bt/b = 1442. 
The horizontal lines indicate the 69 lowest energy eigenvalues, with 30 below 
zero and 39 very close-lying above zero. The inset relates their mean-square 
hyperradii with the absolute values of their energies. The lowest state has 
2mb'^EQ/!f ~ -0.0147 and ^/Wu>/b ~ 136. 



such states are related by 

where the exponential dependence on the strength ^ of the effective potential 
and the number n of the excited state is highlighted. This relation can be 
written as 



En ip') 



!i±i _ „27r/|i| 



En+1 {P^)n 



= e^'^/l^l . (4.12) 



With increasing quantum number these states become exponentially larger with 
exponentially smaller energies approaching zero. 

Around thirty states with this character are obtained for the potential in 
figure 4.5. The lower curve in the inset of figure 4.5 illustrates the relation in 
equation (4.12), and accordingly many states are in this log-log plot represented 
by a point on the straight line with slope —1. The very lowest states deviate 
due to the attraction at small p, and the states close to E = deviate due to 
the external potential. The denser positive energy spectrum in the upper part 
of the inset approaches a straight line with slope -|-1 as expected for a harmonic 
potential. Using equation (2.2) we get 2(p2) = {N - l)(r^2) - 2(iV - l){rf). 
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Even the most bound state with (p^)"^^^ — 1366 then has a root-mean-square 
(rms) distance between two particles (rf2)^^^ — 446, which is much larger than 
the interaction range. Also the rms distance from the centre of the trap (rf )^/^ ~ 
316 is large. 

The intermediate region responsible for the constant A is only present when 
the scattering length is relatively large, i.e. when relation (4.7) is obeyed. This 

corresponds to hypcrradiii larger than p^-^ = N'^^^b and smaller than Pmax = 
N'^^^\as\. The number of Efimov-Hke states A/e located in this region is then by 
equation (4.6) given as 

~ Min (^] 2^ 0.40ArV6in (M) , (4.13) 



where equation (4.9) yielded the last estimate. The number of Efimov-like 
states JVe increases strongly with N. This assumes that the external trap has 
no influence on the hyperradial potential for p < Pmax- However, when the trap 
length 6| is sufficiently small, that is when ptrap — 

extension of the plateau is truncated at large hyperradii. The number of states 
is then estimated by substituting Pmax with ptrap in equation (4.13). This yields 



When the trap length is large and does not terminate the plateau at large 
distances, the mean-square hyperradii of the first and last Efimov-like states are 
of the order p"^-^^ ^ N'^/^b'^ and p"^^^ ^ N'^/^\ag\^, respectively. Equation (4.11) 
then yields the energies of the first and last Efimov-like states 

These energies are independent of the particle number N and remind of the 
kinetic-energy scale of strongly bound two-body states and the two-body binding 
energy, respectively. However, the rms distances f between two particles in these 
many-body states are not given by 6 and a.,. In fact, f contains an additional 
iV-dependent factor, i.e. f ~ A^^/'^6, A^^/"^|as| for the two cases. These constant 
energy limits imply that the density of Efimov-like states increases with the 
particle number. 

These many-body states arise when the two-body scattering length is large. 
This is the condition for the occurrence of the three-body Efimov states [Efi70, 
FJ93], that show characteristic properties similar to equations (4.10), (4.11), and 
(4.12). The author and co-workers therefore proposed to call the many-body 
states with similar attributes for many-body Efimov states [SFJ02a]. A more 
correct name is probably Efimov- /«fce many-body states since some definitions of 
the term "Efimov state" read that infinitely many iV-body bound states occur 
when the N — 1-body system is on the threshold for binding. According to 
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Amado et al. [AG73] "there is no Efimov effect for four or more particles" in 
the sense that being on the threshold for binding in the N — 1-body system 
does not produce infinitely many A/^-body bound states when N > A. This 
statement is not in contradiction with the Efimov-like states discussed here 
since the present Efimov-like A^-body states occur when the two-body system, 
and not the N — 1-body system, is on the threshold for binding. However, the 
quoted remark reminds us that the three-, four-, . . . , — 1-body systems are 
also bound and that many of the A^-body states might be resonances embedded 
in the continua of dimer, trimer, and higher-order cluster states. They could 
be artifacts of the model where only special degrees of freedom are treated, 
and where we recall the possible changes due to larger coupling terms for large 
scattering lengths. However, because the particles are far from each other and 
couplings to the continuum states therefore could bo weak, some of these states 
might be distinguishable structures which could be relatively stable. We return 
to such considerations in section 6.2.4. 

4.5 Trap states and "the condensate" 

In the non-interacting case hyperradial many-body states are located in the 
potential minimum created by a competition between the kinetic energy and 
the external trap. Similar behaviours were seen in the cases of repulsion and 
attraction, see for instance in figure 4.1 the excited states above the lowest two 
located in the minimum at large hyperradii. The corresponding density profile 
of the lowest trap state is similar to that obtained in experiments creating Bose- 
Einstein condensates [DGPS99]. We can therefore call this trap state for "the 
condensate" or define a condensate by the typical signatures of the lowest state 
located in this minimum due to the external trap. 

The attraction produces lower-lying many-body bound states with an av- 
erage distance between the particles much smaller than in the condensate-like 
state. The structure of these states could as well be characterized as a con- 
densate (condensed A''-body state), but they are much more unstable due to 
the much larger density and the larger recombination probability. These lower 
states have no parallels in mean-field computations. 

Through the derived adiabatic potential the two-body unbound mode is 
responsible for the properties of atomic Bose-Einstein condensation where no 
clusterization is allowed. We focus on the state of the condensate in the second 
minimum and in the present work only use the lower-lying negative-energy states 
in connection with the possible decay of the condensate. However, first we 
discuss a definition of a condensate state in the present context. 

4.5.1 A definition of "condensate" 

In mean-field treatments, with repulsive two-body potentials and confining trap 
potentials, the condensate is uniquely defined as a statistical mixture of single- 
particle states with the ground state dominating [PS02, Pou02]. This many- 
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body state is mainly determined by the properties of the trap. It is at best 
only approximately stationary due to the neglected degrees of freedom which 
allow energetically favored di- and tri-atomic cluster states. This instability 
is also an experimental fact seen by permanent loss of trapped atoms, e.g. in 
recombination processes [DCC+01]. 

Without any two-body interaction the properties of the many-body system 
is determined by the thick, dashed potential curve in figure 4.2. Then the 
condensate is a physical state dominated by the ground-state component. With 
attractive interactions (full curve) the deep minimum at small hyperradius is 
produced. Then the ground state, located in this minimum, has nothing to do 
with a condensate. The density is so high that couphngs to other degrees of 
freedom would develop higher-order correlations and processes like three-body 
recombinations would quickly destroy the single-atom nature of the gas. This 
ground state, before or after recombinations, does not show the signature of a 
Bose-Einstein condensate where many particles occupy one single-particle level. 

The formulation in the present work does not use the concept of single- 
particle levels. Therefore we cannot talk about a statistical distribution of 
particles with the majority in the lowest state. However, we can talk about 
a many-particle system described as a superposition of many-body eigenstates 
where the lowest states are favored in thermal equilibrium. To clarify, a quantum 
state is given as the superposition of eigenstates ^n{p, ^) from equation (2.26): 

OC' OC' oo 

^quantum state 

ra=0 n=0 1^=0 

with the normalization X^^ol*^"!^ ~ ^- '^^^^ spatial extension of a conden- 
sate must be sufficiently large in order to exceed a certain minimum interpar- 
ticle distance dc below which the atoms are too close and recombine very fast. 
This distance depends on the scattering length and on the number of particles. 
Therefore, a state cannot be characterized as a condensate if components with 
(rfa) ^ c?c dominating contributions in the wave function. 

One of the stationary states in this model can be defined as the "ideal con- 
densate state", i.e. the state of lowest energy with one component, labeled by 
the quantum numbers Uc and ric, which has 

(r-?2).c,nc > dl . (4.17) 

When no bound two-body states are included in the model, this ideal state is 
determined by the adiabatic component in the lowest angular potential, that 
is fc = 0. On the other hand, the states of lowest energy with v = Q might 
have an average particle distance less than dc. The appropriate choice among 
these excited states depends on the number of particles and on the scattering 
length. The ideal state is then characterized by one dominating component, 
that is z/c = 0, \cnf.\ — 1, and |c„5>£„(,| <C 1. If it is impossible to distinguish 
states with these features, it probably makes little sense to define a condensate. 
The possible states would be too unstable. 
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If d^. is significantly smaller than &t , then the state of lowest energy located 
in the second minimum can be identified as the condensate. This state is charac- 
terized by a radial wave function F{p) with the root-mean-square (rms) radius 
(p^)^/^ approximately equal to the hyperradius at the second minimum of the 
adiabatic potential Uo{p)- 

Figure 4.6a shows the rms interparticle distance f„ given by = (rf2)n = 
2(/9^)„/(iV — 1) for the lowest excited states, labeled by n, in the potential of 
figure 4.2. All states with n < 46 have 206 < r„ < 1006, which implies that the 
particles are well outside the range of the interaction with each other. Whether 
the average distance qualifies a state as a condensate depends on the decay rate 
of this state. From figure 4.6b it is seen that the abrupt change in rms distance 
does not infiuence the energy which changes smoothly for the states in question. 
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Figure 4.6: a) The root-mean-square distance f for = as a function of the 
hyperradial quantum number n for N = 100, as/b = —1, and bt/b = 1442. b) 
The energy E for the same case. 

For the positive-energy states (n > 47) the average particle distance now 
exceeds 20006, that is ~ 36^ which approximately is obtained in the limit 
of a non-interacting gas. This investigation repeated for the higher adiabatic 
potentials ly > 1 results in the same pattern (not shown), although there are 
fewer states with small interparticle distance. 

In section 6.2 we return to a discussion of the appropriate value for dc, which 
then would characterize these states as ideal for a condensate or not. 
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4.5.2 Interaction energy 

The total energy of a state in the first minimum only depends on the interaction 
since this state is bound even in the absence of the external field. Such a 
state has no analogue in mean-field calculations. Total energies of states in the 
second minimum arc dominated by the contribution from the confining field and 
therefore are rather insensitive to anything else than this field. It is then more 
informative to study interaction energies where the large background external- 
field contribution is removed. 

Figure 4.7 shows the interaction energy per particle as a function of the par- 
ticle number for a relatively weak attraction corresponding to the small scat- 
tering length tts/b — —0.84. The two crosses for A'' = 20 (on the left) show 
the results from the two-body correlated model for the lowest adiabatic channel 
u = with quantum numbers n = 7 and n = 8. The interaction energy is nega- 
tive for the lower state, whereas the shown value for the upper state is positive 
due to the extra internal kinetic energy. This is repeated for the larger values 
N = 100 and A'' = 900, i.e. the lowest state characteristic for a condensate is 
shown along with some of the higher-lying states. For N > 950 there are no 
trap states since the barrier has vanished. However, the correlated solutions are 
still stable due to the use of the finite-range potential. 




Figure 4.7: Interaction energy as a function of N for ag/b = —0.84 and bt/b = 
1442. The crosses are results of the present hyperspherical calculation for three 
numbers of particles. The quantum numbers are v = and, for the lowest cross 
in each of the three cases, nc = 7 for iV = 20 {N\as\/bt = 0.012), ric = 52 for 
N = 100 iN\a,\/bt = 0.058), and = 88 for N = 900 (iV|as|/6t = 0.52). The 
sohd curve shows results of the GPE. 

We anticipate the mean-field discussion of chapter 5, i.e. an ansatz for the 
many-body wave function as a product of single-particle amplitudes and a zero- 
range interaction potential lead to the mean-field Gross-Pitaevskii equation 
(GPE). Shown as the solid line in figure 4.7 is the interaction energy obtained 
from the GPE. For small N values the GPE solution is stable and the related 
interaction energy is negative due to the attraction between the particles. A 
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nearly linear behaviour is observed at small particle numbers since each particle 
interacts with N — 1 other particles. We observe the similarity between the 
mean-field energy and the energy for the lowest trap state in the hyperspherical 
model. As N increases, the mean-field attraction increases and a non-physical 
cohapse is inevitable. This instability occurs for A^|as|/6t — 0.55 which cor- 
responds to N 950 with the present set of parameters. There is no stable 
solution to the GPE for N\as\/bt > 0.55. 

Thus, using the correlated model we generally observe both condensate- like 
and collapsed many-boson states. In the present example for relatively few 
particles A^, it was easy to distinguish due to the presence of the intermediate 
potential barrier. In chapter 6 we discuss the case of small or vanishing barrier 
and in chapter 5 compare further to the mean field. 

4.6 Summary 

The radial potential exhibits features from kinetic energy, interactions, and 
external field, and thus combines the information available within the hyper- 
spherical model. The structure of the system depends mainly on the scattering 
length and the trap length. Confined many-body states of negative energy may 
occur even without an external confining potential. This is possible when the 
effective attraction between the bosons is sufficiently large. 

Self-bound states with properties similar to the three-body Efimov states 
occur when the scattering length is very large compared to the range of the 
interaction. These states may still have relatively low density and thus avoid 
the instant collapse due to three-body recombinations. This will be further 
addressed in chapter 6. 

A description of the condensate as effectively a non-self-bound many-body 
state confined by the external trap is possible within this hyperspherical treat- 
ment. Then the average properties are comparable to those of the mean-field 
treatment, as we shall see in more detail in the following chapter. So, the main 
effect of correlations is to allow the Efimov- like many-body states, and the fi- 
nite interaction range prevents an infinite collapse of the system as the density 
increases. 
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Chapter V 



Mean field and validity 



The mean field provides a study of the average properties of a many-boson 
system. Section 5.1 first presents some features of the mean-field method, then 
discusses density-dependent interactions, and finally compares mean-field results 
with the results from the correlated model. Section 5.2 contains a discussion of 
the ranges of validity for both models. 



5.1 Comparison to mean field 

When studying a dilute system of particles, the first approach is usually to apply 
a mean-field method where the many-body wave function for identical bosons 
is factorized into one-particle amplitudes tp as described in section 2.3.1. This 
often leads to the non-hnear Gross-Pitaevskii equation (GPE) [Pit61], where the 
single-particle function enters as a properly normalized single-particle density 
|'0(ri)p. In the stationary case the GPE is written as 



d"^ 1 2 2 47r/i^as,,^ ,,,,, ..o 



V(ri)=0, (5.1) 



where m is the mass of the particles, to is the angular frequency of an external 
trapping potential, and ag is the two-body s-wave scattering length. The eigen- 
value of this equation is the chemical potential /it which is related to the total 
energy E'totai by = dEtotax/dN. 

The GPE was applied to a description of experiments with trapped alkali 
atoms [EB95, BP96] and is now routinely solved for the density profile of a 
condensate [DGPS99, PS02, PSO-S] and for the dynamical evolution of a con- 
densate with the time-dependent version of equation (5.1) [AM02]. Proukakis 
et al. [PBS98] derived a non-linear Schrodinger equation from a microscopic 
treatment of binary interactions and in the low-density limit obtained the GPE. 
The under-lying assumptions arc valid when the scattering length is small 
compared to the inter-particle separations [LS02], i.e. when n|as|'^ <C 1 where 
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n is the density.* The mean-field validity condition is then fulfilled, which 
means that the mean free path ce given by ce = l/(na^) is long compared to 
the average distance which approximately is n~^/^.^ The low-energy scattering 
properties expressed by the scattering length are then decisive. In the following 
we first comment on the choice of interaction, then discuss density-dependent 
interactions, and finally collect some differences between the results from the 
mean-field method and the hyperspherical adiabatic method. 



5.1.1 Two-body interactions 

The origin of the interaction term in the Gross-Pitaevskii equation (GPE) is the 
approximation 

j dr2 V{r2 - ri)|V(r2)|' ^ j dr2 V{ri - rj) (5.2) 

for an interaction of much shorter range than the average distance. This integral 
is given by the Born approximation to the scattering length as we saw in section 
3.1 

/ dr2 V{ri - ra) = , (5.3) 

J m 

where the GPE then occurs with ab replaced by a^. The smallness of a., com- 
pared with the average distance between particles is the criterion of validity 
[Pit61]. This corresponds to a scattering situation where the wave function 
hardly changes due to the scattering, i.e. the wave length is very large corre- 
sponding to low energy and low density. 

The mean-field treatment above corresponds to a zero-range interaction with 
ttB replaced by as 

, , AiTh?as ^, , 

= — ^ — o{r) , r = r2 - ri , (5.4) 

see also section 3.1. This Hmit can be obtained from a finite-range potential 
where the range approaches zero and the strength is appropriately adjusted. 
The finite-range Gaussian interaction of equation (3.7) can be expressed as 

^G(r) = ^^^G(r), (5.5) 
5G{r) ^ -A_e"'-Vb^ , 1= j dr S^ir) . (5.6) 



*In chapters 5 and 6 n denotes the density and not the hyperradial quantum number as in 
chapter 4. 

tThe Danish letter "ce" is pronounced almost like the vowel in the English "them" or 
"Thames". 
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The Gaussian S^^r) is in the Hmit when 6 — > a representation of the Dirac S 
function. For = ob we then have 



However, = Ob is only vaUd when |aB|/& — > 0, which is rarely the case, see 
figure 3.2. 

The aim of computing reliable energies in the mean-field approximation can 
be achieved with equation (5.4) for dilute systems [EG99]. The interaction 
and the Hilbort space must be consistent, that is to say that a renormalized 
interaction follows a restricted space to produce the correct energy. In this 
case the Hilbert space is restricted to the mean-field product wave function. 
Any inclusion of features outside this restricted space, for example two-body 
cluster structures, would be disastrous [FJOla]. Maintaining the finite-range 
interaction with the correct scattering length then results in different properties 
of the interaction even when the range approaches zero on any scale defined 
by the problem. Thus, the mean-field product wave function with a reaHstic 
two-body potential would also lead to wrong results, as we shall see later. 

The full Hilbert space with the correct interaction must produce correct an- 
swers to any proposed question. Whether a realistic interaction combined with 
the present inclusion of two-body correlations reproduces the main features is 
not obvious. However, the investigations in chapter 4 demonstrate that the 
energy in the mean-field approximation for dilute systems is reproduced. This 
asymptotic behaviour is determined by the scattering length which only implic- 
it\y is contained in a given combination of range and strength of the Gaussian 
interaction. This implies that the Hilbert space of the model accounts properly 
for the correlations crucial at large separations. 

5.1.2 Density-dependent interactions 

In section 3.2.4 we studied the angular potential with a zero-range interaction 
as in equation (5.4). This led to the angular eigenvalue Xg from equation (3.31), 
which in the limit of large densities clearly is wrong. A possible treatment at 
large densities is to use a finite-range potential with the correct scattering length 
as discussed in section 3.1. However, this immediately requires a treatment be- 
yond the mean field, e.g. by the Jastrow or Faddeev approaches. Another com- 
mon approach, especially in nuclear physics [SJ87, CBB"'"03], but also for atomic 
many-boson systems [LY57, LHY57, BHM02], is to expand the interaction in 
density-dependent terms. In the present case the two-body interaction can be 
written as 



\imVG{r)=Vs{r) . 



(5.7) 



V{r) = g2{n)6{r) , 



(5.8) 



where the low-density limit of the coupling strength .g2(") is 52(0) = ATiH^as/m. 

We relate the density n to the root-mean-square (rms) hyperradius p = 
a/ (p^). The density is related to the rms separation f, which is defined by 
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= (''12), by 1/n ~ 47rf^/3. The relation = 2p^/{N - 1), obtained from 
equation (2.2), then yields 

3 iV3/2 ^ ^ 

(5.9) 



8\/27r p3 

We then replace p by p and for a fixed hypcrradius calculate the angular po- 
tential as an expectation value of the two-body interactions, equation (5.8), in 
a constant hyperangular function. This yields 

A5„(n)=A5^, A5„(0) = A5. (5.10) 

The reverse relation yields an expression for the density-dependent couphng 
strength, i.e. 

52(n)=52(0)%^. (5.11) 
Aa 

We assume that the numerically obtained angular potentials A, as calculated in 
chapter 3 and parametrized in equation (3.38), represent the density-dependent 
potential rather well, so we identify Xsn with the lowest angular potential for the 
case with no two-body bound states. Here we use the above-mentioned trans- 
lation between n and p. So figure 5.1a shows for various N values \5n{n)/\s = 
g2{n)/g2{0) as a function of the density in the combination N'^n\as\^. 

At low densities the ratio approaches unity which yields the correct limit 
92 — 52(0). At larger densities the deviations are significant, and the couphng 
strength vanishes altogether as n|as|^ — > 00 since the finite-range interaction 
contains no divergence. Figure 5.1b shows the same quantities for different 
scattering lengths and the behaviours are confirmed. 

When the scattering length is negative and large, n|as|^ ^ 1/-^^, the angular 
potential assumes a constant value A ~ Aoo, see equation (3.39) in section 3.4. 
This yields a couphng strength of magnitude 

92{n) Xqc ^, 0.48 ^ 0.77f 

52(0) Xs ~ 7V2/3ni/3|a^| - iV2/3|„^| ' ) 

In this region the couphng strength decreases hnearly with the rms distance f 
between the bosons. The interaction energy per particle in this region is given 
by 

E 1 , 7rft2„2/3 

^-2^2(n)n:.-^^^. (5.13) 

This is independent of the scattering length. A .Jastrow calculation in the denser 
region for positive scattering lengths by Cowell et al. [CHM+02] yielded E/N ~ 
13/i^n^/^/rn, which reminds of the present result for negative scattering length. 
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Figure 5.1: a) The coupling strength in units of the zero-density value as a 
function of the density for ag/b = —10 and various particle numbers, b) The 
coupling strength in units of the zero-density value as a function of the density 
for = 10 and various scattering lengths. The results are obtained by a finite- 
element treatment of the Faddeev-Hke equation (2.75). 



Besides the sign change and the different factor, the present result contains an 
additional dependence on the number of particles. 

This expansion in {n^/'^\as\)~^ is in clear contrast to density-expansions 
where the energy functional is written as expansions in n^/'^jasl [LY57, BHM02]. 
Such a low-density expansion clearly diverges at large densities where it is not 
intended to work. The present results also need corrections at large densities due 
to higher-order correlations, but they might provide a modified zero-range in- 
teraction V{r) = g2{'n)S{r) which could possibly be implemented in a GPE-Hke 
treatment of systems denser than usually within reach. 

In order to avoid collapse due to an attractive two-body 6 interaction, some 
methods apply a repulsive three-body contact interaction. This can be written 
as 

V3(ri2,ri3) = g3{n)S{ri2)S{ri3) , (5.14) 

where we allow a density-dependent coupling strength gs. The expectation 
value of the three-body interactions in a constant angular wave function yields 
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the angular potential 



A3-body(n) = ^N{N - 1){N - 2){N - 3) (iV - ^) (iV - .(5.15) 
Inspired by Gammal et al. [GFTCOO] we parametrize the couphng strength as 

93{n) = ^ , (5.16) 



which for A'' 1 yields, in units of Xs, 

A3-body(n) _ Ntts 



Xs 3[p/{y/Nbt)]^ h 



kz{n) . (5.17) 



Gammal et al. [GFTCOO] use values < fcs < 0.03 which yields X^-hody/Xs < 1 
for a system with p ^ \/iV&t- However, at larger densities we have As.body > ^s- 
Even though this three-body contact interaction can not account for the details 
when three particles approach each other, it might provide a step towards an 
explicit inclusion of three-body correlations."'" 

A treatment of fermion antisymmetry in the hyperangular equations prob- 
ably becomes too complicated when many particles are involved. However, the 
effect of two-body correlations for fermions might be included by a modified 
zero-range coupling strength as described for bosons above. The density de- 
pendence could possibly be extracted for a few fermions and then applied for a 
large number of particles. 



5.1.3 Properties of the wave functions 

In the dilute limit the Hartree wave function is closely related to the hyper- 
radial function and the Jastrow correlated wave function is closely related to 
the Faddccv-likc decomposition of the wave function, see section 2.3. A direct 
comparison of the wave functions is in general not possible as this requires an 
expansion on a complete set of basis functions in one of the coordinate systems. 
The necessary calculations involve non-reducible high-dimensional integrals. 

Instead we use the indirect relations provided in section 2.3.1 where en- 
ergy and average distance between particles are characteristic features of the 
solutions. For a given scattering length the energy E is numerically obtained 
for identical bosons as a function of the particle number. The interaction en- 
ergy is next calculated as E — Eq where Eq = 3Nfiu)/2 is the energy of the 
non-interacting trapped gas. The results for Us/b = —0.84 are shown in fig- 
ure 5.2a. The discussion of stability which follows in section 6.1 shows that in 
terms of variational average distance the GPE energy for attractive potentials 
has a local minimum at large average distances and much lower energies at small 

■tOn this train of thought the work by Bohn et al. [BEG98] was an initiation for the present 
study of two-body correlations. 
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average distances. The physically meaningful mean-field solution is located in 
the minimum at large average distance. This minimum becomes unstable for 
sufficiently large particle numbers. In the example of figure 5.2a no stable mean- 
field solution (solid, thin line) exists for N = 1000. This is consistent with the 
experimentally established stability criterion N\as\/bt < 0.55 [CKT+03] as seen 
from the upper A''|as|/6t-axis. 




Figure 5.2: a) Interaction energy as a function of for a.,/6 = —0.84 {as/b = 
—0.5) and bt/b = 1442. The thin solid line shows GPE results and the plusses are 
obtained from the two-body correlated model; see also figure 4.7. The dashed 
line shows the GPE results for cis/b = —0.5. The upper A^|ns|/6t-axis applies 
for Us/b = —0.84. b) Mean-square distance between the particles for the same 
cases. 



The same figure shows results obtained with the present two-body correlated 
method for three different particle numbers (plusses) . The correlated and mean- 
field interaction energies arc remarkably similar. It may at first appear odd 
that the mean-field interaction energy is marginally lower than by use of the 
correlated wave function which includes an extra degree of freedom. The reason 
is that the mean-field result is obtained with an effective interaction which only 
in the Born approximation has the correct scattering length, while the correlated 
solution is obtained for an interaction with the correct scattering length. The 
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mean-field interaction is effectively more attractive as discussed in section 3.4. 

The proper comparison is then a GPE calculation with ag/b = —0.5, cor- 
responding to the Born approximation for a Gaussian of Ob/6 = —0.5 with 
true scattering length ag/b = —0.84. As seen in figure 5.2a (dashed curve), 
now the mean-field interaction energies are numerically smaller. This compar- 
ison does not include the negative-energy states supported by the attractive 
pocket at short distance, see e.g. figure 4.1. They would appear far below the 
"condensate-like" state shown in figure 5.2a. 

Using equations (2.2) and (2.39), we compare in figure 5.2b (ri2) for the 
solutions to the mean-field approximation and the hypersplierical methods. Also 
this quantity is very similar for the two methods, whereas we again observe 
the discrepancy when we for the GPE method replace Og by Ob. The mean- 
square distance decreases with increasing particle number for calculations with 
an attractive potential. As N approaches 1000, the Gross-Pitaevskii mean-field 
radius approaches zero due to the unavoidable collapse. The same behaviour is 
seen for radii and interaction energies, i.e. the average distance between particles 
decreases until the condensate collapses and the size vanishes in the mean field, 
while many-body bound states with smaller extension play a role in the present 
hyperspherical description. Then also higher-order correlations can be expected 
to be essential and result in recombination processes as will be discussed in 
section 6.2. 

In conclusion, for weak interactions or very small scattering lengths a sta- 
tionary many-body state can be approximated by a product of single-particle 
amplitudes. However, stronger attraction between particles must invoke other 
degrees of freedom like clusterization. Then a single-particle description is not 
valid. This is in agreement with general expectations, and thus confirmed with 
the present point of departure. 

5.2 Validity conditions 

We conclude the chapter by estimating validity criteria for the models. Of 
special importance in relation to trapped atomic gases is a radial wave function 
with rms hyperradius p ~ \/Nbt. Accurate angular eigenvalues in this region 
are therefore crucial for a proper description. If these hyperradii are sufficiently 
large, that \s p <^ VNbt > iV^/^|as|, the angular eigenvalue has reached its 
asymptotic value A ~ Xs- This condition is equivalent to N\as\/bt < N^^^, 
which is obeyed by stable systems with N\as\/bt < 0.55 < N^^^. 

The different models are valid if appropriately designed, i.e. the present two- 
body correlated model reproduces the correct effective interaction for the correct 
scattering length for any short-range interaction, whereas the Gross-Pitaevskii 
equation (GPE) reproduces this same correct effective interaction by using the 
Born approximation. Interaction energies and sizes would be very similar for 
the states corresponding to the condensate. 

From equation (5.9) a given rms hyperradius p is related to the density n of 
the system by n ^ N^/"^ / p^. The zero-range mean-field method (GPE) is usually 
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valid for condensates when nla^P <C 1, see [CHM+02, PS03]. Then the number 
of particles within a scattering volume 47r|asp/3 is on average much smaller 
than one. From the present model, in the zero-range asymptotic region of p > 
A^'^'/^jflsl, we find that np^ > N'^/'^n\as\^ or n\as\^ < l/TV^ < 1, which means 
that the system is very dilute and both the GPE method and the correlated 
method are valid. 

For p < A^^/^|as| the large-distance asymptotic A ~ A5 is not valid. This 
corresponds to that we cannot use a non-correlated model with the zero-range 
interaction. This is here interpreted as an indication of the inadequacies of the 
GPE in the region where p < N'^''^\as\ or equivalently where njagp > 1/N'^. 
This appears different than the usual criterion of validity n|as|^ <C 1, but could 
potentially specify what is meant by "much smaller than unity". 

The present adiabatic liypersphcrical method with two-body correlations 
explicitly allowed in the form of the wave function is valid in the region p > 
N'^/^\as\ where correlations are expected to be insignificant. The inclusion of 
two-body correlations is expected to allow smaller hyperradii p < N'^^^\as\. 
When higher-order clusterizations occur, any method without correlations higher 
than two-body breaks down. The absolute lower criterion must be that the dis- 
tance between two particles on average exceeds the interaction range h, i.e. p > 
N^/'^b. We quote this as the criterion even though explicit calculations might 
prove that higher-order correlations alter the lower limit. 

In conclusion, the validity regions for the two-body correlated method and 



the GPE version of the mean field arc estimated to be 

p > VWb for two-body correlated method , (5.18) 
p > Ar^/^|a^| for GPE (the present result) . (5.19) 

These relations can with equation (5.9) be expressed via the density as 

nb^ < 1 for two-body correlated method , (5.20) 

1 

'"'M < for GPE (the present result) , (5-21) 

n\asf <i; 1 for GPE (usual expectation) , (5.22) 



where we also collected the usual criterion for validity of the GPE. When the 
density is low, both descriptions are valid and the energies are similar. For 
larger densities the importance of correlations increases and the mean-field ap- 
proximation breaks down. At even higher density also two-body correlations 
are insufficient and the particles may clusterize further. 

These conclusions are illustrated in figure 5.3 where the lowest angular eigen- 
value for a case with no two-body bound state and negative scattering length 
is compared to the zero-range angular potential for the same scattering length. 
For low density njogp < the effective energy of the two methods coin- 

cide. For larger densities the GPE energy diverges, while the energy from the 
finite-range model remains finite. Moreover, it deviates in a region where the 
density is still relatively low n|asp < 1 so higher-order correlations, especially 
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Figure 5.3: The lowest angular eigenvalue for ttg/b = —10^, N = 100, and 
no bound two-body states (solid line). The dashed curve is \s for the same 
scattering length. The vertical Hues indicate regions of different density. 

three-body, do not play a role yet. As the average distance becomes smaller, we 
expect corrections due to higher-order correlations. 



Chapter VI 



Macroscopic stability and decay 



Chapter 4 opened a discussion of the macroscopic boson systems known as 
condensates, and a definition of the condensate was given in the present hy- 
perspherical context. In this chapter we discuss stability and consider possible 
dynamics, primarily related to this lowest trap state or "condensate". In terms 
of the degrees of freedom explicitly included in the two-body correlated model, 
section 6.1 presents a criterion for macroscopic stability. This is similar to the 
discussion presented by Bohn et al. [BEG98]. Section 6.2 ventures into a dis- 
cussion of degrees of freedom that in the strictest sense are beyond the present 
model. This especially involves three-body recombination events that poten- 
tially ignite dynamics of the system as one macroscopic whole. 

6.1 Stability criterion 

Macroscopic instability of systems of bosons has been investigated thoroughly 
the last eight years. Macroscopic stability for a Bose-Einstein condensate (BEG) 
means that the BEG state is well defined and has a sufficiently long lifetime 
when considering possible decay modes. In this sense it was originally ex- 
pected that Bose-Einstein condensation could not be reaHzed for atoms with 
effectively attractive interactions. However, this was achieved in 1995 for ''Li 
atoms [BSTH95]. Since then numerous experiments, e.g. [RBG+01, GKT+03], 
have tested the critical region. Also the collapse process itself has been studied 
[SGWH99, DGG+01, RGG+01]. 

The criterion for stability of a system with negative scattering length as 
can be expressed as a critical combination of the number of particles, the scat- 
tering length, and the trap length bt = ^JfiKmuj), where w is a geometric 
mean, sec section 4.1. Recently it was measured that when N\as\/hi > 0.55, 
a condensate of ®^Rb-atoms is unstable [GKT+03]. This can be understood 
as a competition between the kinetic energy which is effectively repulsive and 
the two-body attractive interaction. When the kinetic energy subdues the net 
attraction, a meta-stable system with signatures of a condensate exists. The 
number of kinetic-energy terms equals the number of particles N and the num- 
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ber of interactions equals the number of particle pairs N'^/2. Therefore, as N 
or the scattering length increases, the interactions at some point win and 
for attraction lead to collapse. When the trap length bt is small, the system is 
compressed too much and is according to the criterion unstable. This means 
that the interactions win at larger densities. 

In the mean field this can be formulated variationally with a Gaussian- 
Hartree amplitude with a variational width w, i.e. tp{ri) = exp[— r^/(2i(;^)], see 
also Pethick and Smith [PS02]. The kinetic energy per particle is proportional 
to l/w'^, the external field energy vP /h^, and the interaction energy Nag/w^. 
This leads to the variational total energy 



2m _ 3A^ SA/'w^ 2 



a. 



This yields energy curves in the length scale w, analogous to the hyperradial 
potential curves in chapter 4, which has stable points for a sufficiently weak 
attraction. The critical value is then found to be about 0.67 [DGPS99, PS02]. 
More detailed analysis of the Gross-Pitacvskii equation (5.1), incorporating time 
dependence, yields a value of 0.55 [GFTOl] in agreement with the experimentally 
measured value [CKT+03]. 

By analogy with the mean-field discussion we here take a closer look at the 
derivation of the stability criterion in the hyperspherical frame. This is also 
equivalent to the derivation performed by Bohn et al. [BEG98]. The criterion 
is obtained by estimating when the radial barrier disappears. The effective 
hyperradial potential U{p) from equation (2.29) can in the asymptotic region 
when p > A''^/^|as|, that is when A ~ A5 from equation (3.31), for A'' 1 be 
written as 



If we rescale equation (6.2) with p = ^/3N/2w, this hyperradial potential is 
identical to the mean-field variational energy from equation (6.1), i.e. U{w) = 
Etota]{w). For a large and negative as or a large value of N there is no stable 
region in such a potential. In general, a potential of the form 

u{x) = ^+Bx^-^ , {A,B,C}>0, (6.3) 

diverges to +00 as x — > cxd and if p > 2 to —00 as a; — > 0. For a sufficiently 
small constant C there will always be a local minimum. There will be no local 
minimum when p>2 and 



j5(lJ 2) B(p-2)/4 \p + 2 

For the present case the power p is given by p = 3 and A, B, and C are given 
by comparing equations (6.2) and (6.3). Then the radial potential for negative 
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scattering length has a local minimum only when 



N\as\ 2\/27r 



~ 0.67 . 



(6.5) 



This is identical to the value obtained variationally from the mean field since 
the ag-only dependence at large liypcrradii corresponds to the mean field with 
a zero-range interaction, as discussed in chapter 3 and directly evident by com- 
paring equations (6.1) and (6.2). This gives the right order of magnitude of the 
critical combination of particle number, scattering length, and trap length. The 
discrepancy from the experimental value can be accounted for by the deforma- 
tion of the external field and the zero-point energy from motion in the trap. 
This will be discussed in chapter 7. 

The criterion is in figure 6.1 illustrated by the radial potential, with the an- 
gular eigenvalue obtained from equation (3.38), as a function of the hyperradius 
for a series of different particle numbers and scattering lengths. The strongly- 
varying short-distance dependence is omitted to ahow focus on intermediate 
and large hyperradii. When an intermediate barrier is present, the condensate 
is described as the state of lowest energy located in the minimum at large hy- 
perradius. This minimium exists for ttg < when N\as\/bt < 0.67 as established 
above. 

In figure 6.1a-6.1d the particle number is fixed at iV = 6000 while only the 
scattering length as varies. In figure 6.1d-6.1f the scattering length is fixed at 
tts/b = —0.35 and N is varied. In figure 6.1a the two-body interaction is zero, 
that is Og = 0, which leads to a vanishing lowest angular eigenvalue A = 0. 
The effective radial potential then consists only of the centrifugal barrier and 
the external field with one minimum. In figure 6.1b an attractive potential with 
Us = —0.056 is sufficiently strong to overcompensate for the centrifugal repulsion 
and create a second minimum in the radial potential at smaller hyperradius; the 
final divergence U{p) +oo when p ^ is not included in the scale of the 
figure. An intermediate barrier is left between the two minima at small and 
large hyperradii. A further increase of the attraction in figure 6.1c removes the 
barrier while leaving a smaller flat region. The negative-potential region around 
the minimum at small hyperradius is now even more pronounced. This tendency 
is continued in figure 6. Id with a stronger attraction. With the scattering length 
from figure 6. Id, i.e. a., = —0.355, and a decreasing number of particles the 
intermediate barrier is slowly restored. In figure 6.1e for N = 3000 a barrier is 
about to occur, and in figure 6. If for N = 500 an intermediate barrier is again 
present between a minimum at small and large hyperradii. 

The discussion in this section involved the static properties of states located 
in the local minimum at large hyperradii. Other factors are important when we 
in the following section to some extent include time dependence. 
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Figure 6.1: Radial potentials with bt/b = 1442 and a) iV = 6000, = 0; 
h) N = 6000, as/b = -0.05; c) N = 6000, a,/b = -0.18; d) N = 6000, 
as/b = -0.35; e) iV = 3000, a,/b = -0.35; f)N = 500, a^/b = -0.35. The 
dashed Hnes are obtained with as = 0. 

6.2 Decay 

The Bose-Einstein condensate (BEC) is intrinsically unstable and decays spon- 
taneously, e.g. into lower-lying dimer states. Recombination of two particles 
into a lower-lying state is possible by emission of a photon, but the rate is en- 
hanced when a third particle is involved instead of the photon. This three-body 
recombination process inevitable occurs in a system of bosons when at least one 
two-body bound state exists. This has been suggested to be important for a 
BEC [AdhOl, US03]. The related change of the surrounding medium could lead 
to an instability which involves many particles, and thus result in much faster 
decays which could be described as a collapse [Adh02b]. 

Figure 6.2 illustrates the different behaviours by using the angular eigen- 
values parametrized through equations (3.37) and (3.38). In figure 6.2a the 
scattering length is relatively small and a large barrier separates the outer mi- 
nimum from the inner region. When the scattering length increases, the barrier 
decreases first into a relatively flat region as in flgure 6.2b and then disappears 
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completely as in figure 6.2c when the trap length is exceeded. 
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Figure 6.2: The radial potential from the parametrization for N = 100, bt/b = 
10^, and a) a.,/fe = —6, b) o^s/b = —50, and c) a.,//j —oc'. The shown wave 
function is the lowest radial solution in the non-interacting case. The horizontal 
lines in parts a) and b) indicate an energy level (not to scale). 

The discussion of macroscopic dynamics in this picture involves various iso- 
lated ideas which lead to simple decay rates. These are then incorporated in a 
description of the experimental collapse situation. 



6.2.1 Three-body recombination 

The condensate state is unstable due to the couphngs into degrees of freedom 
different than the coherent many-body mode. The formation of bound-state 
dimers is possible by a three-body process where the third particle ensures con- 
servation of energy and momentum. The number of these three-body recom- 
bination (rec) events per unit volume and time can be estimated by the upper 
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limit given in [NM99, BBHOO]: 

u^ec - 68 , (6.6) 

m 

where n is the density. This expression can be converted into an estimate of the 
recombination rate for a given root-mean-square (rms) hyperradius p = y/ {p^) 
by using the relation between density and rms hyperradius from equation (5.9), 
i.e. n ~ 3N^/'^/{8y/2np^). With the volume V = N/n the total recombination 
rate then is 



■ rec 



N*h\as 



z/,ecV ~ 0.5 . (6.7) 

The recombination rate increases rapidly with decreasing p, as indicated by the 
vertical arrows in figure 6.2 where we interchange p and p to illustrate this effect. 

If N{t) is the number of particles present in the coherent many-body state 
as time t goes by, the recombination time scale Tree can be defined by N{t) = 
N{0) exp(— f/Tree)- This leads to the rate Trec/fi = —dN/dt = N/rrec, so 

_ Nh 2mp^ mf^ 

Since the condensate forms in the external trap, the system must be stable versus 
recombination events on a time scale Ttrap which is given by the oscillator time 
scale 2-!t/u), that is Tree > Ttrap = Stt/w. With the relation l/cj = mb'^/h the 
criterion for a stable condensate becomes 

f > \/8^\a,f/Hy^ = dc . (6.9) 

Here an expression for the minimal separation dc, as introduced in section 4.5.1, 
is obtained. In units of bt we have 

-.^/8^(MY\ (6.10) 



bt \b 



where the determining combination is |as|/6t. 

Thus, for |as|/&t ^ 1 also dc/bt <C 1. The rms distance f for a state located 
in the second minimum is of the order bt and therefore f > dc, i.e. for this state 
f is larger than the critical stability length dc- This state then qualifies as a 
condensate. For *^Rb atoms with Us — 100 a.u. and trapped in a field with 
i^trap — 100 Hz, we obtain Tree ~ 7 days. 



6.2.2 Macroscopic tunneling and recombination 

The second decay process is macroscopic tunneling through the small barrier 
as indicated in figure 6.2b. The model provides stationary eigenstates which by 
definition are time independent. Thus, strictly the states do not tunnel through 
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the barrier. However, an exponentially small tail extends to small hyperradii or 
large density. All particles thus approaching each other would recombine into 
molecular clusters because the density is very large in the inner region. The 
rate of this two-step decay with tunneling through the barrier and subsequent 
recombination is determined by the bottleneck. The rate of recombination due 
to macroscopic tunnehng can be estimated by [BEG98] 



- tunnel 



/■Pout 
J Pin 



dp 



^tunnel — r» 1 



/ 1 d-^Uip) 




m dp^ 


Pmin 



2m 



U{p) - E 



(6.11) 



(6.12) 



where the multiplication by the factor N gives the total number of recombined 
particles. Here Pmin is the position of the second minimum of [/, and and 
Pout are the points where the barrier height equals the energy E. 

When N\as\/ht <C 1, the barrier is large and the very small rate can be 
estimated through equations (6.11) and (6.12). The action integral is then large 
and given by 

cj^lN\n(^\. (6.13) 
2 \N\as\) ^ ' 

Classical turning points of the potential are present when N\as\/ht < 0.53. Close 
to this threshold, i.e. when the barrier is small, the exponent is 

which is valid when A^|as|/(0.535t) is close to unity. 

At the threshold for macroscopic stability then N\as\/bt ~ 0.5, which due 
to the factor of N implies that |os|/6t *C 1. Close to this threshold we have 
f ~ 6t ^ (ic; which means that the average distance between the bosons is so 
large that the three-body recombination is slow compared to the typical time 
scale for oscillation in the harmonic-oscillator trap. Therefore, the three-body 
recombination does not hmit the macroscopic stability of a condensate. In the 
limit (7 <C 1 we get explicitly 

implying that the macroscopic tunnehng process dominates. With cr <C 1 then 
rtunnei/?i - 0.5Ari'tunnei and ftunnei - t'trap, which yields a tuuuehng time of 
about 1/ftrap- For the case with ftrap — 100 Hz the macroscopic tunneling time 
scale is 10 ms. This is much smaller than the three-body recombination time 
scale which close to stability is given by the reciprocal of equation (6.15), that 
is 



~ 7.0N^ > 1 . 



(6-16) 

'tunnel -'- rec 
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The three-body recombination rate is in figure 6.3 shown as a function of 
the hyperradius (soUd curve) and compared with the macroscopic tunneling rate 
(dashed curve) where ah particles in the condensate simultaneously disappear 
during contraction. At small hyperradii the three-body recombination rate is 
much larger than the macroscopic tunnehng rate, whereas the opposite holds 
for large hyperradii. For the chosen set of parameters the two time scales are 
roughly equal around the second minimum where the condensate is located. 
However, the tunneling rate depends strongly on the barrier through the com- 
bination N\as\/bt. Variation of either of the three quantities then moves the 




Figure 6.3: The three-body recombination rate from equation (6.7) in units of 
the oscillator frequency Vtiap = w/(27r), which typicahy is 10-100 Hz, as a func- 
tion of hyperradius for N = 100, a,,/6 = —50, and ht/b = 10*. Shown as the 
horizontal, dashed line is the macroscopic tunneling rate from equation (6.11). 
Shown as the horizontal, dotted line is the macroscopic collapse rate from equa- 
tion (6.17) when the scattering length is much larger than the trap length. 

tunneling rate up or down in figure 6.3. For a larger barrier the condensate 
would only decay slowly by recombination. For a smaller barrier macroscopic 
tunnehng would dominate and ah the bosons in the condensate would partici- 
pate in a "collective recombination" in a short time interval. 

When a few particles rccombinc into dimcrs and leave the condensate, the 
system is no longer in an eigenstate of the corresponding new Hamiltonian. 
An adiabatic adjustment of Hamiltonian and wave function could then take 
place. Since fewer particles and unchanged a., and means a larger barrier, the 
macroscopic stability of the new system is therefore increased. 

6.2.3 Macroscopic collapse 

Scenarios where the boson system develops with time as one unified body are 
open for investigations in experiments where the effective interaction almost 
instantaneously is changed by tuning close to a resonance [IAS"'"98, DCC"'"01, 
RCC"'"01]. An initially small magnitude of the scattering length, corresponding 
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to a stable condensate state in the second minimum of the hyperradial potential, 
can be changed to a value where the barrier is removed. 

In one experiment [DCC"'"01] a condensate is first created with effectively 
zero interaction, i.e. zero scattering length as in figure 6.1a. The radial wave 
function is then located at relatively large distances in the minimum created by 
the compromise between centrifugal barrier and external field. The attractive 
pocket at small distances is not present and the condensate forms as the ground 
state in this potential. Figure 6.4 shows both the radial potential (thick, dashed 
line) and the wave function (thin, dashed line) for schematic model parameters. 

In the experiment the effective interaction was then suddenly changed by 
tuning a Feshbach resonance [CCR+00] to obtain a large and negative scattering 
length [DCC"'"01]. The measurement showed a burst and a remnant of coherent 
atoms. This was interpreted and explained as formation of dimcrs via the two- 
body resonance, a burst of dissociating dimers, and a remnant of an oscillating 
mixture of coherent atoms and coherent molecules [KII02, MSJ02, KGB03]. 

In the present formulation the effective potential is suddenly altered by a 
change of the underlying two-body interaction. The corresponding new radial 
potential, shown as the thick solid line in figure 6.4, has a pronounced attractive 
region which is able to support a number of self-bound many-body states. 




Figure 6.4: Wave functions / and effective hyperradial potentials U in dimen- 
sionless units as a function of hyperradius for TV = 20 and bi/b = 1442. The 
scattering length is zero up to the time t = and then suddenly changes to be 
large and negative at later times t > 0. Potentials and the corresponding wave 
functions are sketched for t = and at a time r ~ 0.1 ms after half a period. 
The horizontal Hues show the stationary negative-energy states for f > 0. 

Since the initial wave function is not a stationary state in the new potential, 
a motion is started towards smaller hyperradii where it would be reflected off 
the "wall" of the centrifugal term. This macroscopic contraction or collapse is 
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indicated by the large arrow in figure 6.2c. If no degrees of freedom beyond the 
model assumptions are involved, the system would then oscillate between the 
centrifugal repulsion and the wall of the external field. This corresponds to an 
oscillation in density. However, during the macroscopic contraction the rate for 
dimer production via the three-body process increases and dimers are produced 
and subsequently ejected from the trap. Since the rate explodes as the con- 
traction culminates, all particles should rccombine instantly. Thus, the relevant 
time scale, i.e. the bottleneck, is the time scale for macroscopic contraction. 

If the only excitations are the degrees of freedom contained in the lowest new 
hypersplierical potential with .s-waves, we can get quantitative information by 
expanding on the new eigenfunctions. The dominating states in this expansion 
are on the transition point between the lowest-lying positive-energy states with 
energies comparable to the initial condensate and the highest-lying Efimov- 
like many-body states, now present because of the large scattering length, see 
section 4.4. These states have a spatial extension as large as that of the initial 
state. The time scale for evolution of the initial state in the new potential is 
then determined by the energy differences between such levels. The states of 
positive energy and large spatial extension confined by the trap are roughly 
separated by the oscillator quantum of energy hco. The corresponding rate for 
populating smaller distances with the consequence of immediate recombination 
is then crudely estimated to be 

^collapse _ 1 ^ j^^^ 



fi Ttrap 27r 

The resulting non-stationary wave function provides a specific oscillation time. 
After half a period the extension of the system has reached its minimum. The 
wave function at this time r rtrap/2 ^ 0.1 ms is also shown in figure 6.4. 
Experimentally [DCC"'"01] the macroscopic-collapse time is verified to be of the 
order ^ 1/uj. These time scales agree on the order-of-magnitude level. 

The rate of macroscopic collapse is also shown in figure 6.3. This is larger 
than the tunnehng rate. The motion in the potential is slow compared to the 
recombination time for distances in the minimum at smah p, whereas the oppo- 
site holds for distances in the minimum at large p. The time evolution after the 
sudden removal of the barrier could then be as follows. A macroscopic collapse 
towards smaller hyperradii sets in. This is followed by emission of dimers which 
lowers the number of remaining particles and results in a reappearing barrier. 
The part of the wave function trapped at large distances in the second mini- 
mum can then stabilize into a condensate with fewer particles. The time scale 
for these processes should then be between the macroscopic-collapse time and 
the recombination time at the second minimum. 

This makes the assumption that no other degrees of freedom are exploited, 
for example the angular dependence of the wave function or molecular bound 
states described by other adiabatic potentials. Direct population of two-body 
bound states requires inclusion of the adiabatic potential asymptotically descri- 
bing these states. This is possible within the model, but constitutes a major 



6.3. Summaxy 



89 



numerical investigation of coherent atoms and molecules, oscillations between 
them, and three-body recombinations within the same framework. Other time 
scales due to these neglected degrees of freedom could possibly turn up in such 
a complete study of the dynamics of a many-boson system. 

6.2.4 Observation of Efimov-like states 

The recombination probability increases with decreasing hyperradius due to 
the higher density, i.e. several particles are close in space and therefore much 
more likely recombine into molecular states. The time scale Tree for three-body 
recombination is given by N{t) = A^(0) cxp(— t/Trec) where N is the number 
of atoms in the condensate. This is as a function of the root-mean-square 
hyperradius p estimated by equation (6.8). This recombination time for the 
highest-lying Efimov-like states with p ~ N'^^^\as\, see equation (3.41), can 
then be compared to the time scale for motion in the condensate which is given 
by Ttrap =2-it/ll). With p ~ N'^^^\as\ we obtain 



Thus, close to the hmit of stability, i.e. N\as\/bt ~ 0.5, we have Tree ^ Ttrap 
for N ^ 1, so the recombination process is rather slow for these highest-lying 
Efimov-like states. Even though the lifetime is shorter than for the initially 
created condensate, it might be long enough for an observation of these states. 

If the Efimov-hke states are populated in experiments where the potential 
suddenly is changed from figure 6.1a to figure 6. Id, they could possibly be 
indirectly observed. A signature of this many-body Efimov effect would be ob- 
servation of the diatomic molecules formed in the recombination process and 
with the estimated time scale from equation (6.18). The rate should then be 
inversely proportional to the square of the scattering length reached after chang- 
ing the potential. The dimers can probably not be distinguished from this and 
other processes, but the measured rate can possibly be separated into different 
characteristic components. Since their lifetime due to recombination processes 
can be very large compared to the time scale defined by the external field, these 
negative-energy self-bound many-body states should essentially maintain their 
spatial extension after the external field is switched off. This is in contrast to 
positive-energy states where only the trap prevents expansion. Thus, a relatively 
slow time evolution of the density distribution without external field should be 
characteristic for these many-body Efimov-like states. A later measurement of 
a system denser than expected for a positive-energy system could then be a 
signature of the self-bound many-body state. 

6.3 Summary 

The stability criterion for the many-boson system is verified within the hyper- 
sperical framework. The avoided collapse at large densities, compared to col- 
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lapse of the Gross-Pitaevskii description, allows a more detailed study of what 
happens during the macroscopic-collapse process. Time scales are estimated 
on the basis of three-body recombination at different densities. Macroscopic- 
tunneling times and oscillation times in a "free-fall" collapse allow the existence 
of the self-bound states without recombination to other cluster-structures. How- 
ever, since the critical region occurs when a zero-range interaction in a non- 
corrclatcd model describes the system rather well, the inclusion of two-body 
correlations does not alter the criterion for macroscopic stability. This is seen 
by considering a trapped gas with n ~ Close to threshold N\as\/bt ~ 1 

implies nlasp ^ This means that n|ns|'' < 1/^^, which is the asymptotic 

region where A ~ and the GPE method is in agreement with the hypersphe- 
rical correlated method, see section 5.2. In this sense the correlations do not 
modify the expectations obtained from a mean-field consideration. However, 
the possibility for a study of couplings between the coherent many-body system 
and the bound two-body channels is clearly beyond the mean field. This is a 
goal for future investigations of decays. 



Chapter VII 



Deformed boson system 



The atom traps in experiments are of cylindrical geometry, as described in sec- 
tion 4.1. For N attractive atoms the stability criterion, as described in section 
6.1, is experimentally established to be N\as\/bt < 0.55 [CKT+03], where 
is the scattering length and bt = \fhj(mw) is the relevant length scale of the 
harmonic trap of geometric average frequency w = ^fuJ^^jT^fol ■ A reduction from 
three dimensions to effectively one or two dimensions was observed experimen- 
tally [GVL+01] in the limit when the interaction energy is small compared to 
the level spacing in the tightly-confining dimension. Experiments with contin- 
uous variation of the trap geometry from three to either one or two effective 
dimensions [GVL+01], with a two-dimensional structure [GBM+01, RENG03], 
and an effective one-dimensional geometry [TOH+03] request a corresponding 
theoretical description. 

Theoretical interpretations and the underlying analyses are frequently based 
on model assumptions of spherical symmetry [BEG98, Adh02b], as discussed in 
section 4.1. Confinement to lower dimensions can also be studied directly with- 
out the three-dimensional starting point. This has been done with a variational 
calculation in Gross-Pitaevskn equation (GPE) [BP96] and more recently in 
the GPE with variational dimensionality [MW02]. Also effects on stability of 
deformed external fields have been investigated by use of the GPE formula- 
tion [BP96, GFTOl, Adh02a,]. Extreme deformations could result in effective 
one-dimensional or two-dimensional systems which can be described by effective 
interactions of corresponding discrete lower dimensions [01s98, PHSOO, PSOl, 
LMDB02]. 

In the present chapter, which corrects and extends the discussions in ref- 
erence [SFJOSb], we rewrite the hyperspherical formulation from chapter 2 to 
account for a general deformation of the external field. Since two-body corre- 
lations are not yet included in the wave function, this hyperspherical approach 
reminds of a mean-field treatment. We investigate the stability criterion in 
section 7.2. Section 7.3 contains an approach to an effective dimension which 
depends on the deformation of the external field. Since the interactions are 
presently not included in this effective dimension, we therefore in section 7.4 
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introduce them on top of the derived rf-dimensional Hamiltonian. Although 
the choice of interactions is not unique, we can with some guess obtain an al- 
ternative stability criterion and subsequently interpret the results in terms of a 
deformation-dependent coupling strength, which is finally compared with known 
results. 

7.1 Hyperspherical description of deformation 

As described in section 4.1 a combination of magnetic fields results in an effective 
trapping potential, which can be described as the deformed harmonic oscillator 
potential Krap acting on all the identical particles of mass m 



The hyperradius p is the principal coordinate, which is separated into the com- 
ponents px, Py, and pz along the different axes, i.e. 



2 2 2,2,2 2,2 2 2,2 n\ 

P = l^Z^rij = Px + Py+ Pz = P±+ Pz , P± = Px + Py^ (7-2) 



where rij = rj—Ti. In the centre-of-mass system the remaining coordinates are 
given as angles collectively denoted by ft, see analogies in chapter 2 and more 
details in appendix A. 3. 

An appHcation here of the method presented in chapter 2 is to assume a 
relative wave function as a sum of two-particle components. In the case of a 
spherical trapping field each two-body component only needs dependence on 
p and the two-body distance Vij = v^psinay through an angle aij. For a 
deformed external field it also needs dependence on the angle between the 
interatomic vector r,j and the axis of the external field. The two-body com- 
ponent should in the cylindrical case then be on the form 4>{p,aij,'dij), which 
would lead to an angular equation in the two variables ai2 and '!?i2 with more 
complicated integrals than those appearing in equations (2.75) and (2.84). We 
will not investigate this, but here restrict ourselves to no dependence on hyper- 
angles. This is expected to dominate for dilute systems where the large distances 
average out directional dependence. 

Thus, we neglect correlations in analogy to a mean-field treatment, so in the 
dilute limit the hyperangular average of the relative Hamiltonian is 



Vtrap(r-i) = -m{u)lx^ + u^yf + ulzf) . 



(7.1) 



N 



i<j 



N 



(7.3) 



2mHq 




(7.4) 
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where d = 1 and b'^ = h/ {mujg) for q = x, y, z. The interactions Vij are averaged 
over all angles fl, which for the zero-range interaction 4Trh^as6{rij)/m from 
equation (3.8) for N ^ 1 yields 



V = 



m 



l<3 



PxPyPz 



(7.5) 



If we replace as Px = Py = Pz = p/V^, this is identical to fi? / {2mp'^) with 
from equation (3.31). 

We define the following dimensionless coordinates and parameters: 



N 
bl + bl 
26? 



_ Px 



b^ 



7 = 



N 

bl 



bz 



26? 



s = 



bxbybz 



(7.6) 
(7.7) 



The deformation along the different axes is then described by (3 and 7, and s is 
the effective interaction strength. The Schrodinger equation HF{px, py, pz) = 
EF{px, Py, Pz) is rewritten with the transformation 



fix,y,z) cx {xyzi''-^^/^F{p„py,p,) 
in order to avoid first derivatives. We then obtain 

1 92 ^ Q2 



(7.8) 



/? + 7 dx"^ /3 — 7 dy'^ 
N'^u{x, y,z) — e 



dz"^ 



2yW^ 



f{x,y,z) = , 



(7.9) 



2 s 
TT xyz 



where e = 2NE/{hjj). Without interaction, i.e. = 0, the ground-state solu- 
tion is 



f{x,y,z) = (xyz) 



{N-2)/1 



exp[-7V(a;2 + + ^^y^^ ^ 



(7.11) 



which for TV ;» 1 is peaked at (x, y, z) = (1. 1, 1). 

Here we do not solve this equation, but instead investigate the character 
of the effective potential u. For axial symmetry around the z axis the x and 
y directions cannot be distinguished, that is when 7 = and [3 = b\/b1 with 
b\ = bxby. This symmetry amounts to replacing p^ and Py by p\/2 in the 
equations. A convenient definition for this case is 2^^ = x^ +y'^. Equipotential 
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b) 




Figure 7.1: Contour plots of u{x,y,z), equation (7.10), with x = y = g as a 
function of {g, z) for s = —0.4/3^/^ corresponding to Nas/b± = —0.4 for three 
deformations. The values for the contours change by 2, 2, and 5, respectively 



for a) /? = 1 (spherical), b) /3 
(pancake-shape or oblate). 



1/16 (cigar-shaped or prolate), and c) /? = 16 



contours of u in the {g,z) plane for g = x = y are shown in figure 7.1 for 
attractive interactions. For a,, < (s < 0) there is always a divergence towards 
— oo when {g, z) (0,0), see equation (7.10). However, a stationary minimum 
is seen in both figures 7.1a (spherical symmetry) and 7.1b (prolate) close to 
{g, z) = (1, 1), whereas this minimum has disappeared for the oblate system in 
figure 7.1c. For very weak attraction a stationary minimum is present for all 
deformations. 

Figure 7.2 shows cuts of the potential u{g, z) along paths close to the bottom 
of the valleys (see inset). The spherical minimum (full line) is shielded by a 
relatively small barrier from the divergence for £i — > 0. The minimum for the 
prolate deformation (dashed curve) is extremely stable although the divergence 
for g ^ Q still exists. For the oblate deformation (dot-dashed line) the local 
minimum has vanished for this attraction strength. 



7.2. Stability criterion for bosons in a deformed trap 
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Figure 7.2: The potential u{x, y, z) for s = —0.4/3^/^ as a function oi g = x = y 
along cuts of the {g, z) plane where z = g^/^ . The height Mbarrier of the 
local maximum at top of the barrier is indicated for the spherical case /? = 1. 
The inset shows corresponding trajectories in the {g, z) plane, compare with 
figure 7.1, for the three deformations. 



7.2 Stability criterion for bosons in a deformed trap 

The barrier height depends on the deformation of the external field, see fig- 
ures 7.1 and 7.2. Extrema {xo,yo,ZQ) of u in equation (7.10) obey the three 
equations obtained from 

|(4-l) = y^^ (7.12) 

K V TT yozo 

and symmetric permutations of x, y, and z. This can be used to determine 
the critical strength s when a local minimum disappears. The results for axial 
symmetry are shown as the thin solid line in figure 7.3. In these units the critical 
strength s is largest for a geometry very close to spherical. Since s = Nag/bt 
and 6^ = bxbybz, this means that at fixed b^, or fixed volume, the scattering 
length can assume the largest negative value for the spherical trap. Gammal 
et al. [GFTOl] performed a time-dependent study with the Gross-Pitaevskii 
equation (GPE) which resulted in the critical strengths here shown as the dotted 
line. This is in large regions lower than the present result, which might be due 
to our neglect of quantum effects and time dependence. The energy gain due to 
correlations is not included in the simple expectation value with the assumed 
angle-independent wave function. The recent value for the experimental stability 
region [CKT+03] is shown as the plus and agrees with the mean-field model. We 
perform two alternative derivations withing the hyperspherical model in order 
to see if the assumptions need modifications. 

These results can be compared to an analytical "spherical" approximation 
where the radial motion is described by only p while the deformed external field 
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Figure 7.3: The critical strength \s\ = A^|as|/6t as a function of the deformation 
/3 = b^/bl from the potential in equation (7.10) (thin solid line), from equa- 
tion (7.16) (dashed line), and from a mean-field Gross-Pitaevskii computation 
by Gammal et al. [GFTOl] (dotted line). The thick solid line is equation (7.18) 
obtained by considering the zero-point energy. The plus is the experimentally 
measured value [CKT+03]. Regions below curves are considered stable in the 
separate treatments. The double- and triple-dashed lines indicate the effective 
cross-overs to two (2D*) and one (ID*) dimensions [GVL+01]. 



^1 
GPE 
deformed 
spherical 
zero-point 
2D* 
ID* 




remains the same. The effective radial potential U is then obtained by adding 
centrifugal barrier and the contributions from zero-range interaction and the 
external field, see equation (2.29). The angular average replaces each of the 
three components and by p^/3 and R^/3, where -R is the centre-of-mass 
coordinates, i.e. 

pj 9 9 

E(^t-p(r.))n = f '-{p' + NR^) , (7.13) 

2mV . A,.. 3/3 



87ra.^(<5(r.,))o = ^yfA^^/^J, (7.14) 



where 51,2^ = h^^ + b^^ + b^^. By comparison with equations (6.2), (6.3), and 
(6.4) the stability condition becomes 

i^<t(ft7), t(l,tl) = ^ = 0.67, (7.16) 
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The spherical Hmit corresponds to 7 = and /? = 1 where the barrier is 
present when \s\ = N\as\/bt < 0.67. The result for a cylinder (only 7 = 0) 
is shown as the dashed line in figure 7.3 and is noticably different from, but 
numerically almost coincides with, the "deformed" treatment, thin solid line. 
An extreme oblate deformation corresponds to the two-dimensional limit where 
bz <C b± and fi ^ 00. Here equation (7.17) yields the critical strength k ~ 
0.4^/0. 61/ 7r/2/?^ As seen from the contour plot in figure 7.1c, the motion is 
now almost confined at 2: = 1. From x = y = g we see that u{q, g, 1) only has a 
local minimum when \s\ < ^/n/2f3~^^^ , which is larger than the value where the 
z motion is not fixed. This is reasonable since more degrees of freedom in the 
model lowers the energy. Baym and Pethick [BP96] obtained with a variational 
study of the GPE the criterion \s\ < ^/n/2f3~^^^ provided that the variational 
width in the axial direction docs not change due to the interactions. This is 
identical to the criterion from studying the potential u{g,g,l), i.e. consistent 
with the fixed value z = 1. This again emphasizes the equivalence between the 
present hyperspherical non-correlated model and the mean-field GPE. 

Analogously, in the extreme prolate limit (one-dimensional) where /3 ^ 0, 
equation (7.17) yields the critical strength k ~ 0.25v^L25\/7r/3^/^. However, 
fixing X = y = 1 in equation (7.10) yields no critical strength since u{l,l,z) 
always has a global minimum. Therefore, the other degrees of freedom are 
essential in this prolate Hmit. 

A modified stability criterion can be obtained by considering the ground- 
state energy Eq of the boson system, which in the non-interacting case is £"0 = 
h{ux + u)y + u)z)iN — l)/2 where the centre-of-mass energy is subtracted. The 
system is unstable when this energy is larger than the barrier height L^barrier, 
the indication in figure 7.2 of the corresponding height Ubarrier for the reduced 
potential. With this condition the criterion of stability is 



where 31^^ = + b~^ + This is seen in figure 7.3 (thick soHd line) to 
be below the GPE calculations [GFTOl]. The improvement is here substantial 
compared to when the zero-point energy is neglected. In particular, for the 
spherical case we get N\as\/bt ~ 0.53 instead of N\as\/bt — 0.67. 

The estimate of equation (7.18) describes the stability problem better since it 
includes the quantum effect due to the zero-point energy. The GPE calculation 
is time dependent and thus describes the dynamics even better and is also 
closest to the experimental value. Since we expect the present non-correlated 
hyperspherical treatment to be in agreement with the mean field, we expect 
that a time-dependent treatment in this frame would yield the same result as 
obtained by Gammal et al. [GFTOl]. 

A recent variational Monte Carlo investigation of the stability criterion in 
elongated, almost one-dimensional, traps yielded the stability criterion niDtfiD ^ 
0.35 [ABGG03b], where tiid ~ N/bz is the density in one dimension and 




(7.18) 
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OiD = —6_l(&j^/os — 1.0326). Equation (7.18) can in the one-dimensional limit be 
written as N\as\/b± < 0.66. The deviation between the two results might be due 
to our use of a three-dimensional zero-range interaction in this non-correlated 
model, whereas Astrakharchik et al. [ABGGOSb] used a one-dimensional model 
with a zero-range interaction with coupHng strength proportional to 1/aiD as 
well as a full 3D correlated model with hard-sphere or finite-range potentials. An 
effective potential analogous to i5(x)/aiD in the general case with intermediate 
deformations would be a rewarding goal. 

According to Gorlitz et al. [GVL+01] the interaction energy is smaller than 
the energy in the tightly-confining dimension when \s\ < a/32/225/3^ for 
the ID Hmit and when \s\ < 1/32/225/3^/3 fQj. ^.j^g 2D Hmit. These cross- 
overs are indicated by double-dashed (two-dimensional) and triple-dashed (one- 
dimensional) lines in figure 7.3. Since the critical region in each limit is below 
the relevant cross-over, stable and strongly deformed systems can be regarded 
as effectively one- or two-dimensional in the sense of these energy relations. 



7.3 Effective dimension 



The deformation of the external field effectively changes the dimension d of the 
space where the particles move. The field changes continuously and d could 
correspondingly vary from three to either two or one. In order to arrive at such 
a description, we aim at an effective d-dimensional Hamiltonian analogous to 
equation (7.4) with only one radial variable p, a deformation-dependent dimen- 
sion d, and an effective trap length bd, i.e. 

2mHa__^_d_,iN-i)-i± ^ ^7 19^ 

fl2 pd{N-i)-iQpP dp^ by ' ^ ' 

where V represents all particle interactions in d dimensions. The requirement is 
that the Schrodinger equation HdGd = E^Gd with rf-dimensional eigenfunction 
Gd and eigenvalue Ed is obeyed, at least on average, i.e. 



/ 



dp p''^''-'^-'G*d{p){Hd - Ed)Gd{p) = . (7.20) 



The lowest free solution, that is with = 0, is given by equation (7.11). 
In the cylindrical case we can relate the rf-dimensional function Gd to this 
by performing the average with respect to the angle 6 in the parametrization 
iP±jPz) = p{sm. 9 , cos 0) . With inclusion of the corresponding volume elements, 
see appendix A. 3, this leads to 

^d(iv-i)-i|(^^(^^|2 ^^3(iv-i)-i r d0cos''-^0sm^''-^0\F{p,9)\^ , (7.21) 

Jo 



where F{p,9) can be obtained by rewriting equations (7.11) and (7.8). 
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The characteristic energy and length can be defined by 

^"=2^,^^-^^^ d62=26i + 62, (7.22) 

which clearly is correct in the three limits, i.e. spherical: c? = 3 and bd = = h±, 
two-dimensional: d = 2 and b±^ b^, and one-dimensional: d= 1 and b^^ b±. 

In general it is not possible to find one p-independent set of constants 
{E(i.bd;d) such that HdGd = EdGd- Instead wc insist on the average condi- 
tion in equation (7.20) with Gd and Ed from equations (7.21) and (7.22). The 
result for axial geometry is a second-degree equation in d with one physically 
meaningful root, see details in appendix E. The results for various N values are 
shown in figure 7.4. The effective dimension depends on N for relatively small 
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/? 

Figure 7.4: The effective dimension d obtained as a function of the deformation 
parameter (3 = b^/b^. Curves for larger N are very close to that for N = 20. 

particle numbers. When N > 20, the curve is essentially fixed. Furthermore, 
the asymptotic values of both d = 1 (small (3) and d — 2 (large (3) are reached 
faster for larger TV since many particles feel the geometric confinement stronger 
than few particles.* Since these effective dimensions are obtained as average 
values over p, the system might look spherical at large distances and strongly 
deformed at small distances, on average resulting in the curves in figure 7.4. 

7.4 Deformation-dependent interactions 

The effective dimension for the non-interacting system possibly changes when 
interactions arc included. The steps of the previous section should in principle 
be repeated with the interactions. However, this would be complicated and miss 
the goal which is a simple effective Hamiltonian with a renormahzed interaction 
in lower dimension, see analogies in the references [01s98, PHSOO, LMDB02]. 

*It may be amusing to speculate on the meaning of d = 2 for elongated cigar-shaped 
confinement (0.1 < /3 < 0.2). 
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We therefore start out with a two-body contact interaction with a coupling 
strength which is modified due to the deformation. This is in line with the 
renormalization in section 5.1.2 due to the inclusion of two-body correlations. 
So we write a rf-dimensional zero-range interaction with a dimension-dependent 
coupHng strength g{d) as 

Vdnj) = 9{d)5^'\nj) , 5(3) = ^ , (7.23) 

where this "d-dimensional 6 function" is defined by d'^'^\r) — for r ^ and 
Jq°° dr r'^~^ S^'^\r) = 1. The distance between two particles, e.g. particle 1 and 
2, is in hyperspherical coordinates defined by ri2 = \/2psina, where the angle 
a enters the angular volume element as 

dn^, = da sin'^-i a cos'^^^^-^^-^ a . (7.24) 

This is vaHd at least for d = 1,2, 3, see appendix A. 2. 3. The effective interaction 
V in equation (7.19) is for ^ 1 then given by the average over all coordinates 
except p: 

y _ Jo^^dn^ VaiV2psma) _ 2N\Nd/AY'^ a. g{d) 

2 2m V{d/2) p^g{3) ' ^ ' ^ 

However, this does not yield instability for c? < 2 since the power d in p~'^ is 
smaller than two, see section 6.1. 

Wc therefore pursue another approach. Inspired by the forms of equa- 
tions (7.14) and (7.25), we write V as 

"^=2^ T{dl2) Jv ' «3 = a«, (7.26) 

which with = coincides with the result for d = 3 if we choose p = 3. The 
effective potential Ud in the rf-dimensional Schrodinger equation corresponding 
to equation (7.19) is then 



2mUd{p) _ 2N^{Nd/A)P/^ ag d^N^ p_ 



2 

(7.27) 

d 



For p < 2 this potential always has a global minimum and thus no collapse is 
present. For p > 2 there is always divergence to — oo when p — > 0. For weak 
attraction, i.e. small \ad\, there is a local minimum. This disappears at larger 
\ad\ when 

jV|a,| ^ hY' 2^+PlH{p-2)(P-^)I^V{d/2) 
h h p(p + 2)(P+2)/4 • ^'-^^^ 

The criterion in equation (7.16) was also obtained by estimating when the crit- 
ical point vanished. Equation (7.16) is valid for all deformations, i.e. any d. In 
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order to be able to compare equations (7.16) and (7.28), we therefore choose 
p > 2 such that equation (7.28) always is applicable. When equations (7.28) 
and (7.16) agree, the efTective interaction strength Orf is given by 

a, _ 2(P-^^/^T{l)5^/\p-2)(P-^y^ , [2TW 

as bt ^/^(p + 2)(p+2)/4/3i/6p/d V 3 ■ ^ ^ 

This effective interaction strength is in figure 7.5 shown as a function of the 
deformation for various choices of the power p. The solid line shows the result 
for p = 3, which is known to be correct for /3 = 1 (d = 3). Similarly the dashed 
line shows the result with p = d, which does not work for d < 2 (P < 0.2). 
Since the effective coupling strength depends strongly on the power p, we need 
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Figure 7.5: The effective interaction strength ad from equation (7.29) obtained 
as a function of the deformation parameter /? = h\/b'j, in the large- limit, 
i.e. the connection between deformation and effective dimension obtained from 
the calculation for = 20 is used for this illustration. The vertical divergence 
of the dashed Hue indicates the inadequacy of the corresponding method when 
d<2. 



further information about how the interactions enter the effective potential. 

An extreme deformation might lead to effectively one-dimensional or two- 
dimensional properties. Pitacvskii and Stringari [PS03] collected results for the 
effective coupling strength in two dimensions that yields g{2) = ^/8Trh'^as/{mbz), 
whereas the result from equation (7.29) in that limit is larger by the factor 
5''/4^4 ^ I Q Even though the results differ by a factor close to two, the right 
combination of lengths shows that we have incorporated the degrees of freedom 
in the correct manner. This was also the case in the previous comparison of the 
stability criterion with the one obtained by Baym and Pctliick [BP96]. However, 
as was also mentioned by Pitaevskii and Stringari [PS03] , in the low-density limit 
in two dimensions the coupHng constant becomes density-dependent, which is 
beyond the present model where correlations are neglected. 

Since p = d for d = I does not yield a meaningful interpretation of the 
stability criterion, a one-dimensional system needs a different treatment. 
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7.5 Discussion 

In conclusion, the hyperspherical method with a non-correlated approach jdelded 
stability criteria as a function of the deformation of the external field. For con- 
stant volume the highest stability was found for spherical traps. Effective di- 
mensions d continuously varying between 1 and 3 were calculated as a function 
of the deformation. The system can possibly be described by a d-dimensional 
effective radial potential with a d-dimensional effective interaction. However, 
this does not have an unambigious form. Applications to restricted geometries 
become simpler, where the obtained two-dimensional coupling strength com- 
pares reasonably with a coupling strength obtained by an axial average of a 
three-dimensional contact interaction. For the one-dimensional case an effective 
coupling strength was not obtained. 

A previous approach to a d = 1 treatment by Gammal et al. [GFTAGO] shows 
that a three-body contact interaction is necessary for the GPE to produce col- 
lapse in one spatial dimension. In the present framework a three-body contact 
interaction for a constant angular wave function produces a hyperradial poten- 
tial proportional to p~^'^, compare with equation (5.15), which for any d > 1 
leads to instability if the thrcc-body coupling strength is sufficiently negative. 
The dotted line in figure 7.5 shows the effective coupling strength for p = 2d, 
corresponding to this three-body zero-range interaction. 

According to Astrakharchik et al. [ABGG03b, ABGGOSa] a Jastrow ansatz 
for a correlated wave function and inclusion of two-body interactions lead to 
collapse in one spatial dimension. According to preHminary Faddeev calcula- 
tions with the two-body correlated model presented in chapter 2, a two-body 
interaction and inclusion of only two-body correlations in one spatial dimen- 
sion do not lead to collapse. It seems that at least three-body correlations or 
three-body interactions are necessary in order to achieve a realistic description 
of collapse in one dimension. 



Chapter VIII 



Conclusions and perspectives 



The present thesis studied few-body correlations within a many-body system, 
especially effects beyond the commonly applied mean field by a method that is 
usually applied to clusterized systems. 

Chapter 2 presented a hyperspherical framework for including two-body cor- 
relations explicitly in the wave function. For bosons this was done as a sum of 
two-body amplitudes, the Faddeev decomposition of the wave function. One 
advantage of this wave function is that it contains a signature of the average 
distance between all particles. In this respect it reminds of the mean field, 
and it is indeed possible to relate the wave functions when the interactions are 
sufficiently weak. On the other hand, in the dilute limit where three-body en- 
counters are rare, this wave function reminds of the .Jastrow factorization into 
two-body amplitudes. Thus, it catches the information from the encounter of 
two particles at the same time as remembering the background cloud of other 
particles. In the shape of a variational equation, the Schrodingcr equation was 
then reduced to a one-dimensional differential equation in a hyperangle, plus 
a simple one-dimensional equation in the hyperradius. A similar result was 
previously obtained by de la Ripelle et al. [dlRFS88] in terms of Faddeev-like 
equations. However, the present equation is variational and the comphcations 
compared to the simpler Faddeev-like equations are not severe. 

In chapter 3 we discussed analytical estimates of the angular eigenvalues. 
These provide results in the dilute hmit in agreement with expectations based 
on mean-field-like assumptions. Then the numerical solutions confirm these re- 
sults in the dilute limit and furthermore provide significant deviations at larger 
densities. In the regime of a very large two-body s-wave scattering length the 
angular potential approaches a constant value which only depends on the num- 
ber of particles. Also the signatures of a two-body bound state are recognized, 
providing a possible link between the scattering channels and the two-body 
bound channels within the many-boson system. 

The macroscopic properties of a trapped system of bosons were investigated 
in chapter 4, where the radial equation was solved for the size scale and total 
energy of the system. Some of the stationary solutions in the hyperradial poten- 
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tial turn out to have a much smaller spatial extension than the typical size scale 
of a Bose-Einstein condensate. Furthermore, the total energy of many states 
might be negative due to a large average attraction at large densities, i.e. when 
the bosons are close to each other. The Bose-Einstein condensate is usuahy 
unstable if the external confinement is removed, whereas these negative-energy 
modes are self-bound, and thus confined even if the trap is turned off. 

Chapter 5 presented the basic assumptions in the mean-field Gross-Pitaevskii 
(GP) treatment of a dilute boson system. We related the obtained angular 
potentials to a possible density-dependent zero-range interaction. This allows a 
possible rcnormalization of the mean-field interaction for calculations with the 
GP equation for denser systems. The energies and average distances for a dilute 
boson system obtained from the mean field are very close to those obtained from 
the two-body correlated model. However, the states with larger densities and 
negative energies have no parallels in the mean-field model. In this region we 
obtained interaction energies independent of the scattering length. Concerning 
validity, chapter 5 gave estimates of the validity ranges, extending the validity 
region for the two-body correlated model to larger densities with deviations from 
the mean field. 

The stability criterion for a Bose-Einstein condensate is in the present hyper- 
spherical treatment of pairwise correlations obtained in a way which is similar 
to the derivation from the mean field. This was evident in chapter 6 where the 
criterion was derived in terms of the liypcrradial potential. Three-body recom- 
bination and macroscopic tunnehng were discussed qualitatively, estimating the 
effects of degrees of freedom that are not explicitly included in the two-body 
correlated model. This is related to discussions by Bolm et al. [BEG98]. The 
discussion of macroscopic collapse when suddenly changing the underlying two- 
body interactions is possible within the model. This provides an estimate of the 
cohapse time which agrees on the order-of-magnitude level with the measured 
time scale. 

A possible improvement is the explicit inclusion of three-body correlations in 
the ansatz for the wave function in order to study three-body recombinations in 
the many-boson system. This can also indicate the validity of the present model, 
and tell if the self-bound many-body states described in chapter 4 have physical 
relevance or if they arc artificial products of the present ansatz with two-body 
correlations. A quantitative study of collapse dynamics can be performed by 
studying the time-dependent problem with inclusion of couplings between the 
different channels in the adiabatic expansion of the wave function. 

The investigation showed that the lifetime of some of the self-bound many- 
body states might be so large that they can be observed in experiments. This 
might be done by turning off the trap, waiting for some time, and then altering 
the two-body interaction such that the particles repel each other and the system 
expands. By extrapolating back to the density profile before expansion, it might 
be concluded that the system did not expand in the time period between the 
external trap was turned off and the two-body interactions were made effectively 
repulsive. 
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The effects of deformation were discussed in chapter 7, which described the 
trapped boson system in a non-correlated hyperspherical frame. The external 
confining field alters the stability criterion in agreement with the experimentally 
measured criterion. For a fixed volume any deviation from the spherical geome- 
try decreases the stability. This was also concluded by Gammal et al. [GFTOl] 
from a study of the time-dependent GP equation. The present treatment pro- 
vides an analytical stability criterion as a function of deformation, which agrees 
well, but not perfectly, with the experimentally measured value. The deviations 
are probably due to the crudeness of the approximations made for the analytical 
estimates, but not due to the lack of correlations since these are also absent in 
a mean-field treatment. We furthermore seeked to use only one length scale to 
describe the deformation. For the non-interacting case this results in an effective 
dimension, which then enters an effective Hamiltonian in a single length scale. 
However, the problem is the inclusion of interactions, which does not appear 
trivial. A proposed effective potential provides a stability criterion, but does so 
in a non-transparent way where the coupling to the two-body interactions has 
vanished. It seems that higher-order correlations or higher-order interactions 
are crucial for a full understanding of stability phenomena in effectively lower 
dimensions. 

An immediate extension of the present work is to complete the treatment of 
two-body correlations in lower dimensions and in the general deformed system. 
This could provide the wanted effective interaction in lower dimensions and 
possibly confirm the results of other approaches, e.g. the references [PHSOO, 
LMDB02] for the two-dimensional cases. 

The treatments in this thesis are performed at zero temperature. The effects 
of a finite temperature can possibly be included as a statistical distribution of 
many-body states, where couplings between the states then will play a larger 
role. This might yield information about the effect of pairwise correlations on 
the condensate fraction and the transition temperature. 

Experiments with trapping of fermionic gases raise many questions about 
the modification of correlations for fermions. Especially the problem of binary 
correlations between identical fermions is a great challenge, but can potentially 
be built on top of a hyperspherical frame. Another approach is an extraction 
of a density-dependent coupling strength for a fcrmion system. This might also 
provide answers to questions in other fields of physics, e.g. in molecular physics 
and nuclear physics, where the mean field is inadequate for studies of exotic 
problems, for example nuclei close to a drip line. 

In conclusion, the present study of two-body correlations yielded insight into 
mechanisms that in the dilute limit can be accounted for by a mean field, and 
yielded deviations especially in the presence of a two-body bound state or a 
resonance. The inclusion of three-body correlations can be the crucial next step 
which provides answers to questions about three-body recombination and the 
structure of lower-dimensional systems. 
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Sammendrag pa dansk (Summary in Danish) 



Resume: Afhandlingen "Korte tilfceldige sammenst0d. Parvise korrelationer 
blandt bosoner" omhandler teoretiske modeUer for parvise pavirkninger mellem 
atomer, der befinder sig i en gas af ens partikler ved meget lav temperatur. 
Teorien bag en model for opdeling af et mangepartikelsystem i sma grupper 
prassenteres, og denne model anvendes i tilfceldet hvor kun to partikler skiller 
sig ud. Specielt unders0ges betydningen af store taitheder og kraftig vekselvirkn- 
ing mellem atomerne. Dette viser afvigelser i forhold til en middelfeltsmodel, 
hvor atomerne ikke har mulighed for at indrette sig efter hinanden. Studiet af 
cigarformede eller pandekageformede systemer indikerer, at man ma inkludere 
pavirkningen mellem tre atomer for at forsta disse systemer til bunds. 

Denne afhandling beskraeftiger sig med korrelationer i bosonsystemer, isaer 
relateret til de mange eksperimenter udf0rt med meget kolde alkaligasser (li- 
thium, natrium, rubidium og caesium) i de seneste 10 ar. Specielt studeres 
afvigelser fra middelfeltet. Efter det introducerende kapitel 1 prffisenterede vi 
i kapitel 2 en hypersfaerisk beskrivelse af fapartikelkorrelationer i et mangepar- 
tikelsystem. Til lavcstc ordcn inkludcrcdc vi topartikclkorrclationcr i form 
af en sum af topartikelamplituder. Nar to partikler kommer taet pa hinan- 
den, ffindrer mangepartikelb0lgefunktionen sig fra den sffidvanlige enkeltpar- 
tikclstruktur, som kcndcs fra cn middelfeltsbeskrivclsc. Pa dcnnc made minder 
den hypersfffiriske b0lgefunktion for et tyndt system om en kombination af en 
Jastrow-beskrivelse og en Hartree-beskrivelse. Schrodinger-ligningen omskrives 
med sadanne topartikelamplituder til en Faddeev-agtig ligning i en hypervinkel, 
der relateres til topartikelafstanden, samt en simpel ligning i hyperradius, der 
beskriver den samlede udstrffikning af mangepartikelsystemet. I den ene hy- 
pervinkel udledte vi ogsa cn mere kompliccrct variationsligning, der scs som 
et alternativ til den Faddeev-agtige ligning. I gr^nsen, hvor middelafstanden 
mellem partiklerne er meget st0rre end den typiske vekselvirkingsrffikkevidde, 
er det muligt at reducere komplikationerne ved denne variationsligning, sa et 
anvendeligt redskab fremkommer. 

I kapitel 3 diskuteredes analytiske egenskaber af vinkelligningen. For et 
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fortyndet system giver dette resultater i overensstemmelse med en ukorreleret 
antagelse for b0lgefunktionen sammen med en vekselvirkning, der normalt an- 
vendes i en middelfeltsbeskrivelse. De numeriske l0sninger bekraefter denne 
grsense, men viser ogsa afvigelser ved st0rre taetheder, bl.a. er tilstedevaerelsen 
af en bundet tilstand mellem to partikler bestemmende for en af egenvaerdierne, 
hvilket ikke forekommer i middelfeltet. Dette kan muligvis danne grundlaget 
for en beskrivelse af koblinger mcllom kondensatfasen og do bundnc tilstande. 

I kapitel 4 sammenf0jede vi effekten fra vekselvirkningerne med signaturerne 
fra en ydre faelde. Dette resulterer i forskellige typer l0sninger til det fulde 
radiale problem. Specielt forekommer mangepartikeltilstande, der er bundet 
selv uden den ydre fseldes indflydelse. 

I kapitel 5 beskrev vi middelfeltsantagelserne, og l0sningerne fra den korrel- 
crcdo metodo sammcnligncdcs mod middolfcltsl0sningcr. Dot or muligt at udlodo 
en tffitliedsafhffingig vekselvirkning fra de korrelerede beregninger. Slutteligt 
vistes det, at gyldigheden af den korrelerede metode strffikker sig til omrader 
med st0rre taetheder og store afvigelser fra middelfeltet. 

Stabilitetskriterier og tidsskalaer for forskellige henfaldsmuligheder diskutere- 
des i kapitel 6, dog uden at vi studerede de tidsafhaengige ligninger, og uden 
at vi inkludorodo koblingor imollom do forskolligo fasor oksplicit. Dotto or on 
mulig udvidelse af metoden. Specielt inklusionen af trepartikelkorrelationer ses 
som en nffirliggende fremtidig unders0gelse. Tilstande med negativ energi kan 
muligvis observeres i eksperimenter, da deres levetid er tilstraekkeligt stor, og 
da deres rumlige udstr^kning vil udvikle sig anderledes i tiden end for systemer 
med positiv energi. 

Deformationons indvirking pa ot bosonsystoms ogonskabor studeredes i kapi- 
tel 7. En ukorreleret fremgangsmade giver analytiske stabilitetskriterier, hvor 
det bedste er i nogenlunde overensstemmelse med et tidsafhffingigt middelfelts- 
studie og med den eksperimentelt malte vaerdi. Vi udledte en effektiv dimension, 
der kan anvendes i studiet af effektive vekselvirkninger for deformerede syste- 
mer eller i en ekstrem graense med diskret, lavere dimension. Formen for den 
effektive vekselvirkning er dog uklar og ma i et n0jere studie af korrelationer 
udledes fra de effektive potentialer. Forel0bige unders0gelser af et endimension- 
alt system viser, at topartikelkorrelationer er utilstrffikkelige for en beskrivelse 
af kendte strukturor. Sandsj^nligvis b0r man inkludere trepartikelkorrelationer 
for at opna en tilfredsstillende beskrivelse. 

Problematikken omkring fermiongasser kan maske udredes med en form for 
topartikelkorrelationer indarbejdet i en hypersfaerisk beskrivelse. Dette kan ogsa 
vise sig frugtbart inden for andre omrader af fysikken, f.eks. i studiet af driplin- 
iekerner eller af molekylaere klynger. 

Umiddelbare udvidelsor af do boskrovno metodor er fa^rdigg0rolse af studiet 
af korrelationer i deformerede systemer og inklusion af eksplicit tidsafh^ngighed. 
Desuden forventes en model indeholdende trepartikelkorrelationer at kunne be- 
svare mange sp0rgsmal, da korrelationer af h0jere orden sandsynligvis er ubety- 
delige selv ved forholdsvist store taetheder. Samtidig kan dette teste gyldigheden 
af resultaterne opnaet vha. antagelsen om topartikelkorrelationer. 
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Coordinate transformations 



From Schaum's [Spi68] p. 124-125 we generalize transformations between large 
sets of coordinates. Wc start with a set of M coordinates denoted by s = 
(si, S2, . . . , sm) and replace these by another set of M coordinates q = {qi,q2, 
. . . ,qM), i.e. Si = Si{q). Wc obtain the relation between volume elements and 
Laplacians as fohows. First define 



dqj 



M 



The Laplacian operators and the volume elements are then connected by 



M 



M 



i=i j=i 

M M 

dsi = H dqj . 



d n d 

dqj h'j dqj 



(A.l) 

(A.2) 
(A.3) 



A.l Jacobi coordinates for N identical particles 

We start with the coordinate vectors for N identical particles in d spatial 
dimensions. These are then transformed to the centre-of-mass coordinates 



1 ^ 



and N — 1 relative Jacobi coordinates ry^ for A; = l,2, ...,A'^— 1: 



N-k 
N-k + 



i( 



N-k 



rN-k+1 



hi:-)- 



N-k 



(A.4) 



(A.5) 
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The inverse relations are 

N-i 



i+l ' 



The notation in relation to equation (A.l) is 

{rix,riy,riz,r2x, ■ ■ ■ ,rNz) for d = 3 , 
{rix,riy,r2x,---,rNy) for d = 2 , (A. 7) 

{rix,r2x,---,rNx) for = 1 , 

and q = {R,r]^_i, . . . ,»7i). The z'th component of the centre-of-mass coordi- 
nates obeys Hr^ = N^/^ and for the relative components /i,,^^ = 1. Then the 
volume element is 

TV N-1 



i=l k=l 

where each vector denotes d degrees of freedom. The total Laplacian is 



i=l 



dr^ tidvl' ^ NOR' 



A. 2 Hyperspherical coordinates 
A. 2.1 Three particles in three dimensions 

For simplicity we first study tho case of hyperspherical coordinates for d = 3 
spatial dimensions and N = 3 particles. The Jacobi vectors are 

(sin cos (/32 \ / sini?i cos<j£'i \ 

sini?2sin(^2 , r]i=pcosa{ sint?isin<pi . (A. 10) 
C0S'!?2 / \ COS-!?! J 

So, we start with six relative coordinates s ~ {rjix, ■ ■ ■ ,i]2z) and wish to obtain 
the volume element and Laplacian operator in the (new) set of hyperspherical 
coordinates q = (p, a, i?2, ¥'2, t^i, <Pi)- In this case hp = 1, ha = p, h^^ — psina, 
h^2 = /3sinasini?2, h^^ = pcosa, and h^^ = p cos a sin i?! , which yields 

7^ p"'' sin^ acos'^ asint?2 sint^i . (A. 11) 

The terms in the Laplacian then become 

'~ ndpi^dp^ p^dp'' dp' ^ > 



A. 2. Hyperspherical coordinates 
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X I d H d 1 d . 2 1 9 /.,„s 

^« = 2^ = 9 . 2 "COS ' (A. 13 

At da da sm a cos^ a da da 

A^, +A^, =- .% , (A.14) 
sin a 

2 15 5 1_ 5^ 

f 

A^,+A^,=- \ , (A.16) 
cos^ a 

A ^ A,o,a, - A« = 4^^'^ - ^ . (A.17) 
Op op 

Al = — — ^ sin^acos^a— + ^^ + — (A. 18) 
sin a cos^ a da da sm a cos-' a 

In this notation hlk is the angular momentum operator associated with ry^. 
A. 2. 2 N particles in three dimensions 

The hyperspherical coordinates are related to the Jacobi coordinates by 

(sini?fcC0S(/3fe \ 
sini?fesinv3fc , fc = 1, 2, . . . , TV - 1 , (A. 19) 
cosi9fe J 

N-l 

pk = Pfc+1 cosa/s+i = p JJ cosaj , p = pjv-i • (A. 20) 

j=fc+i 

The Jacobi coordinates are s = {rjix, ■ ■ ■ ,r]N-i,z)- With ai = n/2 the new 
set of coordinates is q = (p, ajv-i, ajv-25 • • • , 0:2, t^fe, <Pfe)- This yields hp = 1, 
hak = Pk, h^k — Pk sinafc, h^^. = pk sin a/; sini9fc, and the volume element is 

n = p3^-^ • ( n ^^^^M • ( n '^^^^ cos"""^ • 

\ fe=l / V fe=2 / 

Each degree of freedom contributes to the Laplacian as follows: 

"'Hdp hi dp ~ p3Jv-4 0/ dp ' ^^-^^^ 
= ,^jv-i ; — ^sm^afccos^'^-^afe— , A.23 



~ • (nf=fe+i cos2 a,) sin^ a, sin^?, d^^ " dd^ ' ^ ^ ^ 

^ 1 1 1 d^ 
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We then note 

i2 



+ = , ,^N-i ^ ^-4^ ■ (A-26) 



^3 

All the terms can be collected in 



^ = + E + E (^^'^ + J = 4 - %^ > ( A-27) 

i=2 i=l ^ 

A| = n^ + -^ + -4-, kl=i\, (A.28) 
cos^ afc sm 

Afe = - ^sin2 ak cos^'^-^ . ( A.29) 

sm afe cos^'=~* ak ooik oak 



A convenient transformation is 

9A; - 10 



sin ak cos^-''^ "^^/^ ak 

(3A;-4)(3A;-6) , 
+ - T tan afc 



sinafeCos(^'=-^)/^Q!fc . (A.30) 



A. 2. 3 N particles in d dimensions 



Without repeating the steps in the derivation we collect here the results in d 
spatial dimensions. In the general dimension the hyperspherical coordinates 
can be defined in the same way, when apphcable, as for three dimensions. This 
means that for the integer dimensions we have the set of coordinates 

{(p, ajv-i,...,Q!2,</Jiv-i,---,¥'i,'!?jv-i,---,i?i) for d = 3 , 
(p, ajv-i,...,Q:2,<pjv-i,---,<^i) iovd=2, (A.31) 

(p,ajv-i, . • . ,Q!2) ford=l. 

While the coordinates ak G [0, 7r/2] for d = 3 and d = 2, it is for d = 1 convenient 
to include the sign of a Jacobi coordinate in the definition of the corresponding 
hyperangle, and thus the appropriate range for d = 1 is ak S [— 7r/2, 7r/2]. The 
volume element is given by 

dj7fe = dpp'^^^-'^^-^dQN-i , dilk = dQ.^^U^l^J^UQ.k-1 , (A.32) 

fe=i 

dfii = dO^^) , dO(f ) = dak sin'*-^ ak cos'^^'^"^)"^ ak , (A.33) 

fdipkd'dk sint^fc for d = 3 , 
d^k for d = 2 , (A.34) 

1 for d = 1 . 



A. 3. "Hypercylindrical" coordinates 
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The relative Laplacian becomes 

1 — 



A = Ap- 



i 9 d(Af-l)-l ^ 



A2 



pd(JV-l)-l Qp 



k\=l\ 



9 . d-i 
-sin Qfccos 



cos2 ak sin^ ak 

n2 ^ 1 ^ 

sin''"^ afe cos«'('=-i)-i afe Safc 

V«T-sint?fe3^ ford = 3, 

for = 2 , 
for d = 1 . 



(A.35) 
(A.36) 




i'^('=-^^-^afc/-,(A.37) 
dak 



(A.38) 



Useful transformations of the operators are 

= sin(-'-^)/-..cos'^.^^^.. i - ^ + ^'^^'^^^^^'^ + (^-^'^ 
fc^^ cot^ a. + a. cos'-+^ a. , 



+ 



p'<^'"+\ /d,, ^ ^^^^^ . (A.40) 



with non-negative integers ffe = 0, 1, 2, . . , 



A. 3 "Hypercylindrical" coordinates 

Apart from using the same coordinates as in the spherical case, there are, at 
least, two alternative methods for describing a system with cylindrical symmetry 
or different geometries along all three coordinate axes. 



Combination of one and two dimensions 

The relative coordinates, which are important when describing correlations, can 
be described by the usual N — 1 Jacobi vectors, which are now related to two 
hyperradii and corresponding hyperangles by 

(p_L,fesinafc cosy'fe \ 
p_L,fcSinafe sin(^fe 1 , fc = 1, 2, . . . , - 1 , (A.41) 
Pz,k sin Pk J 

where p_L,jv-i = p± and p±^k = P± cosajv-i ■ • • cos ak+i for k = 1,2, . . . , N — 2. 
Analogue relations hold for Pz,k, especially pz,N-i = Pz- The recursions stop at 
/3i = ai = 7r/2. We collectively denote the angles by Q..* The volume element 



*The description can be extended to describe deformations along all axes. Then the px-part 
is separated into px- and py-parts, that are similar to the pz-part in the present description. 
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becomes 



n dn, = dp^p^-^)-^ dpA""-'^-' d^, X (A.42) 
fe=i 

Af-l 

Yl ^dak sin ak cos^^''"^)"^ ak d/3k cos'''"^^"^ dipk ■ 

fe=2 

The relative Laplacian operator becomes 

A = A^, + A^^ ^2 ^ , (A.43) 

where d = 1 for g = = 2 for q =_L, and Ad^jv-i is the operator in d 

spatial dimensions given previously by equation (A. 36). The angles to enter 
equation (A. 36) are for d = 1 the /3fe's, and for d = 2 the a^'s and (fk's ■ 

Parametrization of the hyperradius 

Alternatively, one common hyperradius can be used along with an angle 6 which 
parametrizes the axial and plane contributions as follows: 

sin0 cos QfAT-i • • • cos Q!fc+i sin Qffc coaifk 
Vk = P \ sin6' cosaAT-i • • -cosafe+i sinafe sin(^fe | , (A. 45) 

cos 6 cos Pn-1 ■ ■ ■ cos Pk+i sin Pk 

with Pi = ai = 7r/2. The volume element and the relative Laplacian becomes 

N-l 

J] drik = rfpp3(JV-i)-i^^pos(^-i)-i ^sin2(^-i)-i e X (A.46) 
fe=i 

N-l 

difii Yl [rfafesinafeCos^^''~^^~^afc d/3fcCos^''~^^~Vfe 



k=2 



where A^ jv-i is the operator in d spatial dimensions as before. 
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Hyperangular matrix elements 



B.l Alternative Jacobi trees 



For use in the calculation of matrix elements different Jacobi trees have to be 
chosen [SS77]. The relevant ones in the context of the Faddeev- and angular 
variational equations are shown in figure B.l. 

The coordinates of the standard tree of figure B.la are defined by 



Vn-1 

'nN-2 



Vi = 



N-1 
N 



vn - jY _ ^ (^Jv-i + . . . + n) 



(B.l) 
(B.2) 

(B.3) 



In the (12)(34)-tree of figure B.lb two of the vectors are different from the 
standard tree: 



Vn-2 = ^(^•4-r3), 

Vn-3 = ^{r4 + rs - r2 - ri) 



(B.4) 
(B.5) 



In the (123)(45)-tree of figure B.lc two of the vectors differ from the standard 
tree: 



VN-i 



1. ^ 1. 

-(r5 + r4) - -{ra + r2 + ri) 



(B.6) 



(B.7) 



In the (12)(345)-tree of figure B.ld three of the vectors deviate from the 
standard tree: 

= -^{r4-r3), (B.8) 
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m-i 




Figure B.l: Jacobi trees: a) standard, b) (12)(34), c) (123)(45), d) (12)(345), 
and e) (12) (34) (56). 



r5- i(r4 + r3) 



1/ ^ 1/ 

g (rs + r4 + rs) - -(ra + ri) 



In the (12)(34)(56)-tree of figure B.le four vectors are different: 
»7jv-4 = ^{r4 + r3-r2-r2) , 



'nN-5 



1. ^ 1/ 

-(re + r5) - -(r4 + rs + r2 + ri) 



(B.9) 
(B.IO) 

(B.ll) 
(B.12) 
(B.13) 



Since only inter-relations between TJjv-i, 'nN-2i ^-^d Vn-s needed in 
evaluating the matrix elements, we use the common notation: 



r]N-i = psina , r?jv-2 = pcosasin/3 



(B.14) 



B.2. Matrix elements: Faddeev 
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r?jv-3 = pcosacos/3sin7 , r]^. ■ r]i = rjkili cos 'dk^i ■ (B.15) 

Here i?fe_; is the angle between the /c'th and Z'th Jacobi vectors. Wo abbreviate 
'&N-i,N-2 '&X, '&N~i,N-3 '&y, and '&N-2,N-3 ^ '&z- An azimuthal angle tp 
determining the projection of T7jv_3 o^ito the plane oir]j^_i and t?jv-2 is defined 
in the usual way such that 

cos '&Z = sin '&X sin cos <^ + cos cos . (B . 16) 

With T = {P,'y,'dx,'&y,v} a matrix element of an arbitrary function / of ah the 
variables a and r then becomes 



drfia,r) = l^^^, (B.17) 

~ .cos^^-^S 



/ 

dfg{a,T)= d/3sin2/3cos^^-^°/3 / djsm^ji 
Jo Jo 

X / dt?a;Sint?a; / d^ySm^y / d(p g{a,T) . (B.18) 
Jo Jo Jo 

The normalization is explicitly J dr = 1. In the following matrix elements we 
need relations for intcrparticle distances and therefore define rj^j = {rj — ri)/\/2 
and the angle related to rjij = psina^ = rij/^/2. 

B.2 Matrix elements: Faddeev 

Equations (2.73) and (2.74) are evaluated as fohows. 

In the integral / dr (^<(a34) a convenient choice of coordinates is the alter- 
native Jacobi (12)(34)-tree of figure B.lb. The angle a^A is associated with the 
distance r^A = \/2?734 by the relation 

??34 = W-a = pcosasin/3 = /9sina34 <J=4> sin Q!34 = cos a sin /3 . (B.19) 

The integrand 0(a34) only depends on 034, which is a function of a and (3. 
Therefore at fixed a equation (B.17) reduces to 



/ 



C'^ dp sin2 /3 cos^^-io 13 Hasi) 

dr <^ a34 = ^ 7^ ^ ' = B.20 

sin2/3cos3^-iO/3 

-l=WsI^ / sin2/3cos3^-i°/3</)(a34)^Rir'V(«) • 



To describe three particles in J dr <?i>(ai3) simultaneously, Jacobi vectors of 
the standard tree are needed. The distance between particles 1 and 3 is related 
to the corresponding Jacobi vector 

1 . X 1 ^/3 , , 

Vi3 = -qira-ri) = -VN-1 + —Vn-2 , (B.21) 
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The hyperangle ais, associated with the distance between particles 1 and 3 
through 7713 = ri3/-\/2 = psinaia, is then 

1 3 \/'3 

sin^ ai3 = - sin^ ^ cos^ asin^ /? + sin a cos a sin /3 cos , (B.22) 

where -ffx is the angle between the Jacobi vectors »7jv_i and T7jv_2- Note that 
0(0^13); through q;i3, for fixed a depends on /3 and i?^, which leaves a two- 
dimensional integral. Therefore equation (B.17) becomes 



/ 



dr = ^^^'^ sin- /.cos3^--^o /^/; M,^ sin.. ^(^,3) 

d/3 sin- (3 cos^N- 10 ^ J- sin 

d/3 sin-/3cos3^-i°/3 / dt?. sini?. (/((ais) • 



This integral can be reduced to one dimension by a partial integration. The 
final one-dimensional integral becomes 



/ 



4 r(3jV-6\ 

dr 0(0:13) = /3jv-7\ si^^ ^acos^ ^^ax 

v37rr(^— ) 

/•7r/2-|7r/6-a| 



i 



nTT/ A— |7r/D— a| 

/ dai3 cos^''^"^ 7"'" sinQ;i3 cosai30(ai3) — (B.24) 

^Cq— 7r/3lefa>7r/3l 



(q— 7r/3)e(a>7r/3) 
(7r/3-a)e(7r/3>a) 

rfai3 cos3^~^7~sinai3cosai30(ai3) 



'0 

where sin^ 7=*= = 4(sin- a-l-sin^ Q:i3=Fsinasinai3)/3, and 6 is the truth function. 



B.3 Matrix elements: variational 

We first divide the integrals of equation (2.81) into similar terms, then compute 
them in general, and finally in the short-range limit. 

B.3.1 Numbers of different terms 

We have to evaluate the double sums of equation (2.81) including the potential: 

N N 
k<l i<j 

Three types of terms occur, due to the fact that we vary the wave function 
component 0^2 equation (2.80): the potential concerning particles 1 and 2, 
the potential concerning one of the particles 1 or 2 and a third particle and the 
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potential concerning neither particle 1 nor 2, but a third and a fourth particle. 
We obtain 



N 



N 



N 



N 



^ Vkl = Wl2 + ^ Vil + '^V2l+ ^ Vkl 



k<l 



1=3 



1=3 



3<k<l 



V12 + 2{N - 2)vi3 + -{N-2){N- 3)^34 



(B.26) 



where the arrow indicates the identity of the terms after integration over all 
angles except 0:12, i.e. analogously to the steps leading up to equation (2.72). 
Treating each of these in the quadruple sum, where the repeated use of arrows 
(— >) has the meaning given just above: 

Fixing (j)l2 and V12 yields three different terms: 



N 



N 



N 



N 



i<3 ^ J=3 0=3 3<i<j 



V12 

as shown in figure B.2. 

<l)*,V,(j) 



h2 + 2{N - 2)013 + -{N-2){N- 3)(t>34 



(B.27) 



<l>*,v 



a) 



b) 



3 s 1 



Figure B.2: Illustration of (?i*2^i2-terms. 



Fixing 4>l2 and V13 yields seven different terms. These can be identified in 
two steps, the first of which separates into four different sums: 

N . N N N Af s 

i<j ^i=2 i=3 j=4 4<i<j ' 

Each of these four terms are then identified as: 

JV , JV X 

■^13 Yj ^l-J' ^ '"13 ( + '?^13 + X^ 01j ] 
2=2 ^ i=4 ' 

^ ^^13 {<l>i2 + (/-IS + {N - 3)<^i4] , (B.29) 

VX3 Y = ""13 (4>23^Y ) ^ ^13 [^^'23 ^ - 3)(/)24] , (B.30) 
j=3 ^ i=4 ^ 
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N 

VisY^hj ^VisiN -3)<j>34 , (B.31) 

3=4. 

^ - - 4)<^45 . (B.32) 



fl3 



4<i<j 

The resulting seven types are shown in figure B.3. 




f) ^ " g) 

Figure B.3: Illustration of ^Ij^^is-terms. 

Fixing 4>l2 and V34 yields six different terms, identified as follows. The first 
step is: 

N , N 

■^34 X] = "^34 f ^ 01i + (B.33) 

i<i ^ i=2 

JV JV TV TV 

X] '^2i + X] ^3i + XI '^4i + 51 
j=3 3=4 j=5 5<i<j 

In the next step the sums are treated: 

N . AT s 

V34 X <l>l3 = V34 ( (j)i2 + (pl3 + 014 + X (plj ] 
3=2 ^ 3=5 ^ 
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(B.34) 



V34 4>2j = ^34 ( (t>23 + 4>2A + 02j ) ^ ^^34 [2(/)l3 + {N - 4)</)i5] ,(B.35) 
N 



VSi X] 



j=4 

N 



N 

V34 \(I)34 + Y^ hj ] [034 + {N ~ 4)(^35] , (B.36) 

j=5 



^34 X] '^4i ^ W34(^ - 4)(/)35 , 

^ 1 
^^34 51 -^^y ^ f^34 ^(^ - 4)(iV - 5)(;i.56 • 
5<i<j 

See the six types in figure B.4. 



(B.37) 
(B.38) 



a) 



c) 



e) 



b) 



d) 



f) 



Figure B.4: Illustration of 022^34-terms. 



B.3.2 Evaluation of terms 

Tfie term of figure B.2a is trivial since the integrand is independent of r. The 
terms of figures B.2b B.2c, B.3a, B.3b, B.4a, and B.4c can be evaluated by 
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equations (B.21) and (B.24). 

The term of figure B.3c becomes with the use of the standard Jacobi tree of 
figure B.la 

J dr /(ais) 3(023) = -^f^m^g^ ^ 

/■7r/2 pTT 

/ sin2/?cos^^-i°/3 / rfi? sini? /(aia) 5(a23) , (B.39) 

Jo Jo 

sin^ q;i3.23 = — cos^ a sin^ ^ 4 ^^'^^ '^^^ ^'^^ ^ ^^'^ ^ '^'^^ ^ ' (-^ -^0) 

The term of figure B.4f becomes with the use of the alternative (12)(34)(56)- 
tree of figure B.le 

2^jv 

dr v{azi) <l){a56) = x (B.41) 

TT 

j dp sin2/3cos3^-^°/3 j d-y sin'^ ^cos^^''^^ via^i) ^(ase) , 

sin a34 = cos a sin /3 , sin ase = cos a cos /? sin 7 , (B .42) 

^jv = (3iV- 8)(37V- 10)(37V- 12) . (B.43) 

The terms of figures B.3g, B.4b, and B.4d are evaluated using the (123) (45)- 
and (12)(345)-trees of figures B.lc and B.ld, so 



J dT h{a,T) = ^ d(3 sin2/3cos3^-^°/3 X (B.44) 

/ d'j sin^7COS^^"^^7 / d'&x,z sin ^J^^^^ Is (a, r) , 
Jo Jo 

where /5(a, r) can be either v{a34)^{a35) or /(ai3)s'(a45). The relevant angles 
are 

sin = cos a sin /? , sin 045 = cos a cos /? sin 7 , (B .45) 

. 2 COS^ a / 2/1-2 -2/1 

sm a35 = — ^ — ^^3 cos p sm 7 + sm p 

+2\/3cos/3sin/3sin7cosi?2^ , (B.46) 

and q;i3 given by equation (B.22). Note the identity J dr v{as4)(j){ai5) = 
Jdr ^(ai3) v{a45)- 

The terms of figures B.3d, B.3e, B.3f, and B.4e are evaluated using the 
standard Jacobi tree. Then equation (B.17) reduces to, with i = 1,2,3, 

drfiais) g{a^i) = ' dp sin^ /Jcos^^-^o d-d^ sini?^ x 

dj sin2 7cos^^-"7 / d'dy sin'dy / d(p /{aia) gia^) ■ (B.47) 
Jo Jo 



I 
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The angles aij can be determined by p sin aij = rjij through the relations 

\/3 1 

r 

f2 1 1 



Via — ~^Vn-2 + 2^N-\ ) 



1 



V2i 



1 



(B.48) 
(B.49) 

(B.50) 

(B.51) 



B.3.3 Results in the short-range limit 

The integrals in the short-range limit, when the range h of V{rij) is much smaller 
than the size scale p, are: 



/ 
/ 



/ 
/ 
/ 
/ 



/ 



dr V 
dr V 
dr V 
dr V 
dr V 
dr V 
dr V 
dr V 
dr V 
dr V 



a34)<A(Q;56) - vi(a)R^^"^V(Q!) , 
Q:i3)'/>(ai3) = Ri^~^^?;0(a) ^ V2ia)(j){0) 
ai3)^{ai4) V2{a)R\£4(f){a) ; i = 1,3 , 
Q;i3)'/>(a23) - V2{a)(p{a) , 

Q!l3)<^(a24) ^ t;2(Q!)Ri324 > 

Q;i3)(/)(a45) - ?^2(q;)R! 



The integrals are given by 

4 r(3i^) 



3jV-12^ 



7r/2 



d7sin 7 cos' 



3iV-13 



7 <^(«2i) 



(B.52) 
(B.53) 
(B.54) 
(B.55) 
(B.56) 
(B.57) 
(B.58) 
(B.59) 
(B.60) 
(B.61) 



(B.62) 
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where sin 0:35 = \/3 cos a sin 7/2, sin = cos a cos /Jo sin 7, sin/3o = ta,na/\/3, 
and 



2 T{^) 



Jo 



d7sin 7 cos 



3JV-13 , 



sin^ a?4 = - sin^ a + ^ cos^ a cos^ /?o sin^ 7 



2V2 



H ^ sin a cos a cos /?o sin 7 cos i?x , 

3v3 

sin^ = ^ sin^ o; + ^ cos^ a cos^ /Jo sin^ 7 

H p sin a cos a cos /Jq sin 7 cos ??x , 

3v3 

sin^ a?o = - sin^ a + - cos^ a sin^ 7 
"4 2 

+ — p sinacosasin7cos??x . 



(B.63) 



(B.64) 



(B.65) 



(B.66) 



" (2) 

The two-dimensional integral ^ijl-i^Pioi) can be reduced to a one-dimensional 
integral, analogously to equation (B.24), by a transformation of the general 
form 

/■7r/2 /•TT 

1=1 d'ysm^^cos^^ I Msm.^ (j){a') , (B.67) 
Ja Jo 

sin^ a' = f{a) + g'^{a) sin^7 2f{a)g{a) sin 7 cos 1? , (B.68) 
to the one-dimensional integral 

r '+ 

1 r 

' ~ ' ^ 'MAXjO.a'-} 

da' sin 2a' cos^"*"^ 7" ^(0;') 



2f{a)g{a){p+l) 

f.MAX{0,-a'-} 



(ia'sin2a'cosP+^ 7+ (/.(a') 



. ± ±sina'-/(a) . ,± ti \^ 1 \ 
sm7 = — , sma = j{a) ± g{a) . 



9{ol) 



(B.69) 
(B.70) 



The function MAX{a;, y} outputs the largest of the two numbers x and y. 
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Properties of Jacobi functions 



The Jacobi function V^'^\x) is related to the hyperangular kinetic-energy 
eigenfunctions. More specifically, an eigenfunction to the operator Il\ from 
equation (A. 37) is vi°'''^\x) with x = cos2Q!fc, a = (d — 2)/2, and b = d{k — 
l)/2 — 1. Some relevant properties of these functions in the relation to the hy- 
perspherical treatments are given by Nielsen et al. [NFJGOl]. The important 
properties in this context are the following [AS65, NFJGOl]. 

The Jacobi function vl""'''\x) is a solution to the differential equation 

(1 + x'^)y"{x) + [a + b+{a + b+ l)x]y'{x) +u{u + a + b+ l)y{x) = . (C.l) 

A second solution is T',?'"' (— .r), which for integer v is identical to vi°"'''\x). 
For non- integer ly the Jacobi function 'Pv{x) is regular at a; = 1 and irregular at 
X = —1. We will for integer v not consider the irregular solution, which diverges 
at both X = ±1. 

Important relations for the Jacobi function are 

^^'•°'<'" ° r(. + i)r(-. _ + - + + 1: 1 + °; - 

V{a)V{v + b+l) (1 + 



T{v + l)T{-v)T{nu + a + b+l) 



xT(^ + b+l,-v-a;l-a; , (C.3) 

where T is the hypergeometric function 2^1 [AS65] . 

For small values of the argument ak they can be rewritten via 

Vi^'^\x) = Vl'''''\-x)cosTrv-Q';^'''\-x)smTru , (C.4) 
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Q(-'')(cos2a) = r(a)r(. + 6+l) ^^^^^^ 
" ^ ' itT{v + a + h + l) ^ ' 

Some properties of the gamma function r(x) are r(a; + 1) = xT{x), 
r(a;)r(l -x) = 7r/sin(7ra;), and T{A + x)/T{A) for A > \x\. 

Rotation properties 

For integer values of the quantum number u the Jacobi function for d = 3 and 
fc = — 1, i.e. a = 1/2 and b = 3N/2 — 4, can be written as the polynomial 
(omitting upper indices) 

P^(cos2a) = ^c^,„sin2"a, (C.7) 

n=0 

c.,0 = 1, (C.8) 
(3jV-5)(^-»)+2(^^-n^) 
c.,n+i - c,,„ (n+i)(2n + 3) ' ^^'^^ 

Some rotation properties of these kinetic-energy eigenfunctions are 

R(^-2)po = Ri4"'^Po = n, (C.IO) 

R(^-2)p, = -^^1- (C.12) 
Related properties are 

R^r'^sin^a = -^_L-^[3+(7V-5)sin2a] , (C.13) 

Rir^sin-a = ^;I^-^^|±^cos-a. (C.14) 



Appendix D 



Numerical scalings of angular potential 



Chapter 3 contained an account of the properties of the numerically obtained an- 
gular eigenvalues. We collect here some of the details behind the parametrization 
of the angular eigenvalue in equations (3.37) and (3.38). They were published 
as a part of a larger article [SFJ03a] and are kept here for completeness. 

D.l Effective dependence on the scattering length 

The angular eigenvalue spectrum coincides with the free spectrum at both small 
and large hyperradii; at p = because all interactions are multipHed by and 
at /5 = oo because the short-range interaction has no effect at infinitely large 
distances. Thus, perturbation theory for a Gaussian potential shows that for 
small p the eigenvalues all change from their hyperspherical values Ai/(0) = 
2v{2v + 3Ar - 5) with z/ = 0, 1, ... as 

A.(p) - A.(0) = ^^^iV(7V - l)p2 . (D.l) 

If the two-body potential is attractive, but too weak to support a bound 
state, the eigenvalues reach a minimum and then return to one of the finite hy- 
perspherical values. For a two-body bound state of energy one eigenvalue 
diverges as A = 2mE^'^^ / fi? . The corresponding structure describes, appropri- 
ately symmetrized, one pair of particles in that bound state and all others far 
apart from the pair and from each other. In addition to this finite number of 
such eigenvalues the hyperspherical spectrum emerges at large distances. 

To illustrate we show in figure D.l a number of angular eigenvalues A as 
functions of hyperradius for different potentials. The entirely positive (solid line) 
corresponds to a repulsive Gaussian. The curves diverging at large hyperradii 
(dotted and thick dot-dashed lines) correspond to potentials with one bound 
two-body state. 

The convergence of A as p ^ is due to the finite range of the potential and 
depends on the interaction range h. The deep minima at small to intermediate 
distances depend strongly on both the number of particles and the strength of 
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10 



1 1 1 
10* 



102 



10^ 



10^ 



p/b 



Figure D.l: Angular eigenvalues X,^ (divided by 10^) as functions of hyperradius 

divided by interaction range, p/b, for N = 100, for scattering lengths a,,/6 and 
numbers of bound two-body states A/b indicated as \v{as/b,N'B,) on the figure. 



the attraction. Increasing the strength of the attraction leads to larger negative 
values of A. This trend continues by increasing the attraction even further until 
the same scattering length is reached but now with one bound two-body state. 

The asymptotic behaviour of A is compared to the zero-range result A5 in 
figure D.2. The convergence to the limiting value is fastest for the smallest 
value of as (dashed curve) already reflecting that the correlations arising for 
large scattering lengths (dotted line) cannot be accounted for by the simple 
zero-range result. This is well understood for three particles where the Efimov 
effect (very large as) extends correlations in hyperradius to distances around 
four times the average scattering length [FJ93, JGF97]. These effects are not 
present in the mean-field type of zero-range expectation value contained in A5. 
When p exceeds as by a sufficiently large amount Xg is approached. 

The positive scattering length also leads to an eigenvalue approaching A^ 
at large distance with a similar convergence rate (soHd curve). A stronger at- 
traction corresponding to one bound two-body state produces one diverging 
eigenvalue while the second eigenvalue converges towards the lowest hypersphe- 
rical value (dot-dashed curve) . It almost coincides with the lowest eigenvalue for 
the same scattering length but for a potential without bound two-body states 
(dotted curve). 

The deviations from Xs at large distance is in all cases less than 10%. The 
asymptotic Ijehaviour is very smooth but still originating in systematic numer- 
ical inaccuracies which can be cured by increasing the number of integration 
points in the finite-difference scheme. 
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Figure D.2: Same as figure D.l, but the angular potential is shown in units of 
the zero-range result in equation (3.31). 



D.2 Dependence on the number of particles 

The angular eigenvalues increase rapidly with N as seen from the A^^/^-de- 
pendence in A^, equation (3.31). The major variation in magnitude is then 
accounted for by using this large-distance zero-range result as the scaling unit. 
Figure D.3 shows a series of calculations for the same two-body interaction for 
different numbers of atoms. Ah curves are similar, i.e. there is a systematic 
increase in the characteristic hyperradius, where the curves bend over and ap- 
proach the zero-range result. The large-distance asymptote is determined by 
the scattering length. A characteristic length pa is conveniently defined by 

Pa{N) = TV^/'^la.l , (D.2) 

where the power is obtained numerically to be very close to the indicated value 
7/6. 

The quality of this scaling is seen in figure D.4 where all curves essentially 
coincide for distances smaher than pa- At larger hyperradii the zero-range re- 
sult of -hi should be obtained. However, here numerical inaccuracies produce 
systematic deviations from a common curve, i.e. the deviations increase with 
N. 

The numerical curves can be rather well reproduced by the function 

A(-)(iV,p) = \\s{N,p)\-g(-Hp/pa) , (D.3) 
g(-)(x)=,geo(l-e-/-)(l + ^) , (D.4) 

where ^oo has the value —1 in accurate calculations, because < 0. The 
exponential term reproduces the rather steep approach to the asymptotic value 
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Figure D.3: The lowest angular eigenvalue as a function of hyperradius for 
a,/b = -401 for four different numbers of particles N = 10^, 10^, 10^, 10^ The 
angular potentials are in units of . 

as seen in figure D.4. The behaviour at smaller distance, depending on the 
range of the interaction, is here simulated by the Xb-term. The extreme limit 
of /o ^ is not computed and not included in the approximate function in 
equation (D.4). 

The two groups of computations in figure D.4 are reasonably reproduced by 
the parameter sets Xa — 0.74, Xb ~ 2.3 • 10~^, and goo — —0.8 or g^o — —0.4 for 
the high and low accuracy, respectively. These parameters may also depend on 
the scattering length. Table D.l gives the best choice of parameters for different 
as- 



tts/b 


-5.98 


-401 


-799 


-4212 


-Qoo 


0.99 


0.80 


0.65 


0.30 




1.06 


0.74 


0.59 


0.28 




0.93 


1.081 


1.099 


1.077 


Xb 


0.15 


2.3 • 10-^ 


1.15- 10"^ 


2.2 • lO-'' 


Xb/{b/\as\) 


0.92 


0.922 


0.919 


0.927 



Table D.l: Numerical values of (/ooi Xa, and Xb for four scattering lengths. 

It should be noticed that — ,9cc and Xa both are of order unity, and that 
the fraction Qoo/xa is almost constant, except for the smallest scattering length. 
The parameter xt, introduced to account for the finite interaction range, is close 
to h/\as\. 

At large hyperradii, where Xa p/ pa or equivalently p ^ N'^^^\as\, A*^"' 
approaches Qoo | A^ | . The expected large-distance asymptotic behaviour is A*^~ ^ 
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Figure D.4: The same as figure D.3, but with p in units of pa- The larger points 
following the intermediate curve {goo = —0.8) are calulated with the highest 
numerical accuracy and the smaller points along the upper curve (goo = —0.4) 
are obtained with lower accuracy. The curve for g^o = —1.0 is the expected 
correct asymptotic behaviour. 



Xg and goo should therefore approach —1 in increasingly accurate calculations. 
The results for N = 100 and different scattering lengths, see figure D.2, confirm 
this conclusion by deviating less than 10% from Xs at large hyperradii. 

A well-established result for iV = 3 identical bosons is the large-distance 
behaviour [JGF97] 

A,(iV = 3,p) = i^, {D.5) 

which is in agreement with A^- obtained from equation (3.31) for = 3. Then 
the universal function g^~^ asymptotically approaches goo = —1 for all scattering 
lengths. 

The function g*^"^ is almost independent of This combined with the 
conclusion for iV = 3 impHes that ^oo = — 1 is vaHd for all scattering lengths 
and particle numbers. 

The angular eigenvalue is given by g^^^{x) ^ gocx/xa for Xb ^ p/ Pa ^ Xa. 
Numerical calculations in this intermediate region of hyperradii therefore rather 
accurately determines the fraction goo/xa — —1.08 as given in table D.l. With 
goo = —1 this implies that Xa — 1/1.08 ~ 0.92. The parameters of g^~\x) in 
equation (D.4) can now be collected to be 

5<x> = -l, a;„~0.92, a;6 ~ 0.92 . (D.6) 

\as\ 

The accuracy of the parametrization is seen in figures D.5a-d, where the 
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angular eigenvalues are shown in units of A^") with the individual set of pa- 
rameters from table D.l. A fairly good agreement is found for p/pa > xi,. 
The remaining deviations occur at small hyperradii, which is not included in 
the ^("^-parametrization, and at large hyperradii where the numerical inac- 
curacy increases with increasing scattering lengths. On the other hand the 
large-distance behaviour is known from analytic considerations, which renders 
numerical computations at these distances superfluous. 




o 



1 

0.8 
0.6 
0.4 
0.2 





IQ- 



TV = 10^ 
TV = 10^ 
N = 10^ 
= 10^ 

* □ X * 

J , l_lJ , 



(d) 



_L 



10- 



IQ- 



10-1 

P/Pa 



10° 



IQi 



10^ 



Figure D.5: The lowest angular eigenvalue A in units of A^"-*, equations (D.3) 
and (D.4) and table D.l, as functions of the hyperradius in units of pa, equa- 
tion (D.2). The scattering lengths are given by a) Us/b = —401, b) tts/b = —799, 
c) as/b = —4212, and d) Us/b = —5.98. The different A''-values are as indicated. 
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D.3 Bound two-body state 

In the presence of a bound two-body state one angular eigenvalue eventually 
diverges at large hyperradii as 

A(2)(p) = , i;(2)<o, (D.7) 

where E^^^ is the energy of the two-body bound state. In the Hmit of weak 
binding, or for numerically large scattering lengths, the energy of the two-body 
bound state is given by 

E(^) = -^c, (D.8) 
ma% 

where c approaches unity for large scattering lengths. 

The angular eigenvalue corresponding to a two-body bound state is para- 
metrized by an expression similar to equations (D.3) and (D.4). The effect of 
the bound two-body state only shows up at large distances where the behaviour 
corresponds to equation (D.8). The small and intermediate distances resemble 
the behaviour when no bound state is present. Therefore the angular eigenvalue 
is given by the parametrization 

X^+HN,p) = \Xs{N,p)\g^+\p/pa) , (D.9) 

»'^'(^)-('-?)fe-4\/f«0' <°-'°' 

with the notation and estimates from equation (D.6), i.e. Xb — 0.926/|as| and 
gcc /:i'a — —1.08. The terms in the second bracket of this expression only aim at 
the correct behaviour in the limits of small to intermediate and large hyperradii. 
The exact transition between these regions is not reproduced. 

Figure D.6 shows a comparison of the parametrization in equations (D.9) 
and (D.IO) with the computed angular eigenvalues for a potential with one 
bound two-body state. For the large scattering length in figure D.6a one smooth 
curve applies for all the particle numbers; numerical inaccuracies set in at larger 
hyperradii, which is most obvious for the largest particle numbers. This smooth 
curve is in a large interval of hyperradii at most deviating by 20% from the 
parametrized form, and even less than 10% at large hyperradii before numerical 
instabilities set in. 

The shape at intermediate distances could be improved for example by inclu- 
sion of a linear term in equation (D.IO). The smooth curve at small hyperradii is 
outside the range of validity of the parametrization, i.e. this is within the range 
of the two-body potential and therefore depends on details of the interaction. 

The lowest eigenvalue diverges at large hyperradius as described in connec- 
tion with figure D.6a and equation (D.9). If the two-body potential only has 
one bound state the second eigenvalue is expected to approach zero at large 
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Figure D.6: a) The lowest angular eigenvalue Aq in units of A^^^ , equations (D.9) 
and (D.IO), for Us/b = +100 and c 1.02, when the potential holds one bound 
two-body state. The number of particles is indicated on the figure. The pa- 
rameters are Qoo/xa = —1.09 and Xi, = 9.2 • 10"*^. b) The first excited angular 
eigenvalue Ai in units of Aj for ag/h = +10. 



distances as Xs- This pattern should be repeated with more than one bound 
two-body state, i.e. the first non-divergent angular eigenvalue should behave 
as A^ for large p. Figure D.6b therefore compares the computed first excited 
angular eigenvalue with A^ for different N. As in figure D.4 smooth and al- 
most universal curves are obtained at small p, where the approach to unity 
sets in exponentially fast depending on N, but now for p one or two orders 
of magnitude larger than pa . A parametrization would also here be possible. 
The large-distance asymptotic behaviour of the first excited state corresponds 
to a repulsive potential as seen by the approach to -|-1. However, at small and 
intermediate hyperradii the potential is attractive (Ai < 0). 



Appendix E 



Derivation of effective dimension 



For the case of N non-interacting identical bosons trapped in a cylindrically 
deformed iiarmonic field the Hamiltonian can be written with the choice of 
coordinates in appendix A. 3 as 

H = Hp + He + fn, (E.l) 
_ 1 d 3jv_4 d , p'^ 

~ cos^-2esin2^-3^^5e de 

Here Th is the angular kinetic energy operator which is neglected later when 
there is no dependency on the internal angles. The ground state total wave 
function is known to be 

N / 2 



with hq = yjhj (m.ujq). We change to the coordinates {p. 0} given by pz = pcos 
and p± = psinO, see appendix A. 3, and obtain the wave function 

NRl NR]_ 

^2 



*totai = exp ( - ^ - — 5^ ) F{p, 6) , (E.5) 



F{p, 9) = exp ( - ^ ) exp 



„2 „„„2 , 



1 1 



(E.6) 



We write the Schrodinger equation and integrate over the coordinate 6 as follows: 
= j deQe{e)F*{p,e)(Hp + He-E)F{p,e) , (E.7) 
OeW = cos^-^sin^^-^^e. (E.8) 
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Doing this we end up with terms only depending on p. Inckiding formally the 
integration over hyperradius, but without completing it, we get 



= Jdp p3(Jv-i)-i J ne{e)F*{p, 9) (ilp + He- E^F{p, 6) 



dp -E+—{N- l)(2Ci + a) 

L Im 

, r(a)r(2a) cn-'.J ^ . 2\ N-l 



(E.9) 
(E.IO) 



where AA{a,b,z) is the Kummer function, i.e. identical to the confluent hyper- 
geometric function \T\{a, b, z) [AS65], Cq = 1/6^, and Ai = — C_l. Further- 
more, it is clear that 



E = ^{N - l){2Ci_ + C) = noj^{N - 1) + ]^huj,{N - 1) 
We desire to write an effective d-dimensional Hamiltonian as 



2m 



d 



pd{N-i)-i Qpi- 



d(N-l) 



dp bi 



(E.ll) 



(E.12) 



with an effective dimension d, and a general length scale bd- The normalization 
of the corresponding d-dimensional Schrodinger equation, with eigenvalue E^., 
is 

J dp p''(^-i)-iGS(p) - Ed) Ga{p) = , (E.13) 

where Gd is a d-dimensional wave function. We want to approximate the correct 
equation (E.9) by this equation. From the two normahzations we identify the 
d-dimensional wave function Gd by the relation 



p3(^-i)-ii?(p)=p'^(^-i)-i|G,(p)|^ 
Using this in equation (E.13) we obtain 

dp p^^-^'^-^Gliip) [Hd - Ed) Gdip) 



(E.14) 



(E.15) 



/ 



dp fi 



,3(iV-l)-lj 



/l2 



2m 



v{x) 



1 



-x{Cz ~'^d)--i^- d)a[a{S + d) - 2] 
-^C^A^xpi-A^x) + ^Ajxp^i-A^x) , 



f \ _ M{a,'ia + 1, z) 
" M{a,3a,z) 



(E.16) 
(E.17) 
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If we subtract equation (E.15) from equation (E.9) we get 

J dp [E-Ea + ^v{p')] = . (E.18) 

Wc aro now in a position to study the three hmits i) spherical {bz = b±), ii) 
two-dimensional {b± ^ bz), iii) one-dimensional {bz ^ b±). 

When the external trap is spherically symmetric, i.e. bz = b±, we have 
Cz = C±, Ai = 0, and /i(0) = 1. This yields 

vip') = -p\Cl - Cl) - ^(3 - d)[a(3 + d)-2]. (E.19) 

The bracket of equation (E.18) is zero for all hyperradii when three conditions 
are true: 

Ed = E=^{N -l){2C^+Cz) = \hw{N -I) , (E.20) 
Cd = Cz ^ bd = bz = b^ , d = 3 . (E.21) 

E is the ground state energy minus the centre-of-mass energy for identical 
non-interacting particles of mass to in a three-dimensional oscillator of frequency 
w. 

The two-dimensional geometry occurs when the external trap is squeezed 
along the r-plane such that b± ^ bz- This leads to Cz ^ C± and Ai = 
Cz — C± ~ > 0. In this case the hyperradius is determined by the radial 
trap length, i.e. p ^ ^/Nhj_, which implies that typically — Aip^ ~ —NCzb\ = 
—NCz/C±_ —oo. We therefore need the limit of ii{z) when z ^ — oo: 



a 2a} — a 

z 

This leads to 



3 / a 2a? — a\ ,^ ^ 

,{z):.-(l^- + ^^). (E.22) 



v{p') ~ -2aCz - p\Ci - Ci) - -^{2 - d)[a{2 + d)-2]. (E.23) 
The bracket of equation (E.18) becomes to order 

^ " '^'^ + " ^"^^ ~ ~ ^^^^ ~ 7^'^ ~ "^^t"^^ + - 2]} • (E.24) 

This is zero for all hyperradii when 

= E- ^2aCz = - 1)2<^± =nw±{N-l), (E.25) 

Cd = C^ ^ bd^b^ , d = 2. (E.26) 

The system is effectively one-dimensional when the external trap is squeezed 
along the ^;-axis such that bz ^ b±. This leads to C± ^ Cz and Ai = Cz — C± ~ 
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— C_L < 0. In this case the hyperradius is determined by the axial trap length, 
i.e. p ~ VNb;,, which implies that typically -Aip^ ^ NC±bl = NC^/C^ -» 
+00. We therefore need the limit of n{z) when z — > +00: 

K.).f(l-^i^). (E.7) 

This leads to 

v{p^) ~ -4aCi - p\Cl ~ Cl) - -^(1 - d)[a{\ + d) - 2] . (E.28) 

The bracket of equation (E.18) becomes to order 

E-Ea + ^[-AaC^- p\Cl - Cl) - ^(1 - d)[a{l + d) - 2]} . (E.29) 

This is zero for all hyperradii when 

Ed = E~^AaC^ = ^{N-l)C,= '^nuj,{N-l), (E.30) 
Cd = ^ bd = h, d=l. (E.31) 

The results are collected in table E.l. In all three cases the energy is given 



limit 


spherical 


oblate 


prolate 


condition 


bz - b± 


b± > bz 


bz > bjL 


z = -Aip' 





—00 


+00 


p{z) 


1 




3q. I'-j^ g— 1 \ 


Ed/{N-l) 






\tkJz 


bd 


bz ^ b_i_ 


b± 


bz 


d 


3 


2 


1 



Tabic E.l: The (typical) values of z, fi{z), Ed, bd, and d in the three limits: 
spherical, oblate, and prolate. 

by 

Ed = ^fuVdiN - 1) , cOd = ^. (E.32) 
2 mb^ 

In the general case, when we cannot assume bz = b± or bz ^ b± or bz <^b±, 
it is not possible to obtain a simple expansion in p of the bracket in equa- 
tion (E.18). Instead we study the integrated equation. 

We rescale the equation in a convenient length scale bo given by 



bl EE 2bl + bl 



(E.33) 
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In the three special hmits we found that Ed = h'^ad/{mb'^) and 6^ = bl/d. We 
therefore introduce the parameters and Pd given by 



Ed 



h?aed 



(EM) 



This leads from equation (E.18) to the integrated equation 
1 



= erf--/3i+5(rf)/i(/3)-/2(/3) 



= a'l[x-'] , (E.35) 



MP) = 2Pr +Pz- \PI + - Pr)I{x^lJ^[x\Pr " /?.)]} 

+ ^(/3.-/3r)'/{xV[^'(/3r -/?.)]}, 

g{d) = (3-(i)(3 + d-2/a) , 
bl 



P 



bl 



P =^-2+1 

Pr- ^2^-'^+ 



nf{x)] = / dx x'^-' exp{-Prx') 

a,3a,x''{pr-Pz)]f{x) . 



xM 

If we use the expectations that 

ed = d^l+e), l3d = d{P + 5)^'^ 
we get that 



£ < 1 



(5 < 1 



Ad^+Bd + C = 0, A=l/2 + e-5/2- h{p) , 
B = 2h{p)/a , C = (9 - 6/a)MP) - f2{P) • 

If we demand only one solution, i.e. B"^ — AAC = 0, we get 

B 2C 
^ = -2A=-^- 



(E.36) 
(E.37) 

(E.38) 

(E.39) 
(E.40) 



(E.41) 
(E.42) 



(E.43) 



The results for various iV-values are shown in figure 7.4. 

Furthermore, we check that £ — 6/2 is smah (< 10^-^) so the approximations 
in equation (E.40) are valid. The conclusion is that generally we can use Ed = 
h'^ad/{mb'^) and 6^ = bg/d as the relevant energy and length scales, respectively, 
with d given by figure 7.4. 
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Appendix F 
List of notations 



IVnf at inn 


T^pcpri nf inn 




tts 


two-body s-wave scattering length 


3-7 




Born approximation to Ug 


2,3,5 


ad 


d-dimensional interaction parameter 


7 


am 


one-dimensional scattering length 


7 


b 


interaction range 


2-6 


bt 


trap length 


2-7 


bd,bq,bx,by,bz,b± 


trap lengths 


4,7 


d 


spatial dimension 


7 


dc 


condensate length 


4,6 


E, En 


total relative energy 


2,4-7 


^'total 


total energy 


5,6 




two-body energy 


3 


Ed 


energy in d dimensions 


7 




reduced hyperradial wave function 


2,4,7 


foo 


Efimov hyperradial wave function 


4 




hyperradial wave function 


2,4,7 


9,92,93 


coupling strengths 


5,7 


G 


kernel for variational angular equation 


2 


Gd 


d-dimensional wave function 


7 




reduced angular Hamiltonian operator 


2,3 


h 


Planck's constant 


2-7 


H,Hq 


Hamiltonian operator 


2,7 


Hd 


Hamiltonian operator in d dimensions 


7 


i 


complex number i = \/— 1 


2 


i, j, i 


indices 


2-7 




Boltzmann's constant 


4 


K = Kn-1 


grand hyperangular momentum 


2-4 


Kk 


hyperangular momentum 


2-3 


In.k 


generalized angular momentum 


4 
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Notation 


Description 


Chapter 


h 


angular momentum 


2 


It 


thermal length 


4 


li, I2 


trap lengths 


7 


Lk 


collection of angular momenta 


2 


L = ijv— 1 


total relative angular momentum 


2 


£ 


Laguerre polynomial 


4 


m 


particle mass 


2-7 


TOfc 


projection of angular momentum 


2 


M 


total mass 


2 


Mk 


projection of 


2 


M = Mjv-i 


projection of L 


2 


n 


single-particle density 


2 


n 


density 


3,5,6 


n 


hyperradial quantum number 


4 


N 


number of particles 


2-7 


M 


number of bound states 


4 


A/b 


number of bound two-body states 


3 




number of Efimov-Hke states 


4 


6 


operator 


2 


Pi 


momentum 


2 


P 


permutation operator 


2 


Pb 


total momentum 


2 


V,V 


Jacobi function 


2,3 


q 


degree of freedom (e.g. x, v, z) 


4,7 




coupling term 


2 


r = ri2,rij 


interparticle coordinates 


2-7 




particle coordinates 


2-7 


r, fn 


root-mean-square (rms) distance 


4-6 




rms separation from mass centre 


4 


R, Rq 


centre-of-mass coordinates 


2,7 




rotation operators 


2,3 


ReS 


effective range 


3 


S 


reduced interaction strength 


7 


t 


time 


6 


T 


temperature 


4 


f 


kinetic-energy operator 


2 


Tp 


hyperradial kinetic-energy operator 


2 


u 


reduced two-body wave function 


3 


u 


reduced hyperradial potential 


6,7 




hyperradial potential 


2-4,6,7 


V, Vij 


reduced two-body potential 


2,3 


Vl,V2 


rotations of potential 


2 




two-body potential 


2-3,5,7 
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Notation 


Description 


Chapter 


V 


total interaction potential 


7 


Vd 


two-body potential in d dimensions 


7 


^rap 


external trapping potential 


2,4,7 




Gaussian two-body potential 


3,5 


Vi 


zero-range interaction potential 


3,5,7 




three-body interaction potential 


5 


V 


volume 


6 


w 


variational width 


6 


x,y, z 


cartesian axes 


2-7 


Y 


spherical harmonics 


2 


y,y 


hvperspherical harmonics 


2 


ce 


mean-free path 


5 


a = ajv-i = ai2, ex' 


primary hyperangle 


2,3,7 


Oiij 


hyperangle related to Vij 


2,3,7 


Oik 


hyperangle related to ?7fe 


2,3 


/3 


deformation parameter 


7 


7 


deformation parameter 


7 


r 


Euler's gamma function 


2,3,7 




decay widths 


6 


5 


Dirac delta function 


3,5,7 




d-dimensional delta function 


7 




Gaussian representation of 5 


5 


A 


Laplacian operator 


A 


£ 


reduced energy 


7 


f]k 


Jacobi coordinates 


2 


9 


polar angle for hyperradius 


7 


e 


step or truth function 


2,3 




polar angle 


2,7 


K 


wave vector or wave number 


2,3 


A, Aj/ 


hyperangular eigenvalue 


2-6 


Ak 


kinetic-energy hyperangular eigenvalue 


2-4 


Xs 


A with expectation value of Vg 


3-7 


Acx) 


plateau value for angular eigenvalue 


3-6 


A(2) 


angular potential for two-body state 


3 


As-body 


A from three-body interaction 


5 


AjV-l 


grand hyperangular momentum operator 


2 


h 


hyperangular momentum operator 


2 




chemical potential 


5 


= VN-UVk 


hyperangular quantum number 


2-4 


^rec 


three-body recombination rate 


6 


^tunnel 


macroscopic tunneling rate 


6 




trap frequencies 


3-7 




Efimov scale 


4 
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Notation 


Description 


Chapter 




partial hyperangular momenta 


2 


p 


hyperradius 


2-7 


On-, P \ , Pr 


components of hyperradius 


7 


p 


root-mean-square hyperradius 


3-6 


Q 


planar projection of coordinate 


7 


a 


action integral 


6 


dr 


reduced hyperangular volume element 


2,3 


T T T 1 r 


time scales 


6 


T 


centre-of-mass wave function 


2 




Faddeev component 


2,3,7 




Faddeev kinetic-energy eigenfunction 


3 


4> 


reduced Faddeev component 


3 




angular wave function 


2-4 


$K 


angular kinetic-energy eigenfunction 


3 




azimuthal angle 


2 


X 


phase shift of two-body wave function 


3 




wave function component 


2,5,6 




total wave function 


2,4 




angular trap frequencies 


2-7 


0, i^TV— 1 


all hyperangles 


2-4,7 




part of hyperangles 


2 




volume element for ak 


2,7 




angles for direction of rik 


2 
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